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ABSTRACT 


Boundary  layer  convective  eddies  (BLCE)  observed  via  synthetic  aperture  radar 
(SAR)  on  17  June  1993  during  the  Hi-Res  2  field  experiment  are  investigated  using  a 
three-dimensional  nonlinear  spectral  model  derived  from  the  shallow  Boussinesq 
equations.  Boundary  conditions  are  appropriate  for  marine  atmospheric  boundary  layer 
(MABL)  flow  and  allow  for  heat  and  momentum  fluxes  across  the  lower  boundary,  as  is 
required  to  determine  sea-surface  stress  patterns.  Enhancements  made  to  the  model  for 
this  study  include  the  addition  of  an  inversion  layer,  an  unstable  surface  layer  and  weighted 
nonlinear  numerical  dissipation.  A  more  numerically  accurate  integration  scheme  is  also 
incorporated. 

A  linear  stability  analysis  of  Hi-Res  2  data  finds  maximum  growth  rates  for 
disturbances  having  horizontal  wavelengths  that  match  those  determined  from  the  17  June 
1993  SAR  image.  Desirable  by-products  of  the  analysis  include  the  determination  of 
physically  plausible  mean  horizontal  wind  profiles  and  inversion  strengths. 

Two  model  integrations  are  performed  using  identical  measured  atmospheric  input 
values  and  initial  conditions,  but  with  different  dissipation  weights.  The  run  with  weaker 
dissipation.  Case  1,  exhibits  a  fully  three-dimensional  cellular  structure  that  is  consistent 
with  SAR  imagery  and  that  is  aligned  nearly  normal  to  the  direction  of  the  mean  wind. 

Case  2,  with  stronger  dissipation,  is  quasi-two-dimensional  and  aligned  14®  to  the  right  of 
the  direction  of  the  mean  wind.  Cross  sections  of  potential  temperature,  streamfunction 
and  vertical  velocity  for  both  cases  reveal  that  the  inversion  is  capping  the  BLCE  with 
only  minor  anomalies  in  the  inversion  layer.  Magnitudes  of  several  dimensional  quantities 
are  compared  with  observations  and/or  other  studies  with  encouraging  results. 
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Chapter  1 
INTRODUCTION 

Since  its  civilian  introduction  nearly  three  decades  ago,  satellite-  and  aircraft-bome 
S5nithetic  aperture  radar  (SAR)  has  become  a  viable  instrument  to  aid  the  study  of  the 
marine  atmospheric  boundary  layer  (MABL).  The  potential  usefulness  of  SAR  to 
meteorology  and  oceanography  was  predicted  by  Vesecky  and  Stewart  (1982),  and  SAR 
today  continues  to  show  great  promise  for  MABL  research,  as  recent  studies  suggest 
(e.g.,  Alpers  and  Briimmer  1994;  Sikora  et  al.  1995;  Mourad  1996). 

SAR  works  by  detecting  variations  in  radar  reflectance  from  the  ocean  surface 
(Thompson,  et  al.  1983).  The  backscatter  patterns  are  generated  by  differences  in  surface 
roughness  and  are  mapped  to  create  high-resolution  images.  Wind-imposed  stress  is  the 
primary  mechanism  for  the  generation  of  centimeter-scale  capillary  waves  (V esecky  and 
Stewart  1982),  and  this  wavelength  magnitude  corresponds  to  that  which  reflects  waves 
via  Bragg  resonance  scattering.  Generally,  SAR  images  depict  high-stress  regions  as 
white  and  low-stress  regions  as  black  (see  Sikora  etal.  1995,  for  a  high  quality  SAR 
image).  In  the  absence  of  intense  oceanographic  and  synoptic-scale  atmospheric  wave 
forcing  phenomena,  we  may  assume  that  the  waves  and  resultant  sea  surface  stress 
patterns  are  due  almost  entirely  to  smaller  kilometer-scale  atmospheric  effects,  namely 
secondary  boundary  layer  circulations  (Sikora  et  al.  1995).  Therefore,  SAR  appears  well 
suited  to  studying  the  wind  variability  caused  by  secondary  circulations  within  the  MABL. 
By  proper  interpretation  of  a  SAR  image,  we  can  theoretically  deduce  the  conditions  in 
the  MABL  that  create  the  observed  stress  patterns,  as  well  as  regions  of  convergence  and 
divergence  at  the  surface  (Gerling  1986;  Alpers  and  Brummer  1994).  Before  employing 
SAR  as  a  reliable  remote  sensing  tool  for  the  determination  of  surface  layer  wind  speed 
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End  direction,  we  must  first  fiilly  understEnd  the  effects  of  MABL  flows  on  the  seE  surfsce 
(Thompson  et  al.  1983).  To  Eccomplish  this  objective,  here  we  model  secondEry 
circulEtions  in  the  MABL  End  study  their  effects  on  modulEting  seE  surfEce  stress. 

There  Ere  two  mEjor  types  of  secondEry  boundsry  lEyer  circulEtions  thEt  CEn  develop 
End  orgEiiize  the  seE  surfEce  stress  pEttems.  The  first  comprises  quEsi-two-dimensionEl 
feEtures,  commonly  referred  to  es  boundEry  lEyer  rolls;  the  second  comprises  three- 
dimensionEl  flows,  CElled  boundEry  lEyer  convective  cells  (Brown  1980).  Three- 
dimensionEl  cells  mEy  form  independently  of  rolls  when  sur&ce  lEyer  winds  Ere  light,  or 
fi'om  rolls  when  thermEl  forcing  is  increEsed  sufficiently  (Woodcock  1940, 1975; 
DeErdorff  1976).  There  Ere  three  principEl  instEbility  mechEnisms  responsible  for  the 
formEtion  of  rolls  End  cells  in  the  Etmospherei  inflection  point  instEbility,  pETEllel 
instEbility  End  thermEl  instEbility  (Brown,  1980;  Stensrud  End  Shirer,  1988;  Etling  End 
Brown,  1993).  The  first  two  drEw  energy  from  the  bEckground  wind  sheEr  in  the 
boundEry  kyer  (pEller  1965).  The  inflection  point  instEbility  mechsnism  tEps  the  normEl 
component  of  the  meEn  sheEr  End  SEtisfies  the  REyleigh  End  Fjortofl  conditions  by 
requiring  e  point  of  inflection  End  e  vorticity  mEximum  within  the  bEckground  wind  profile 
(Kundu  1990).  PsTEllel  instEbility  drEws  energy  from  the  pErsllel  component  of  the  menn 
sheET  viE  the  Coriolis  terms  (Lilly  1966).  Due  to  the  length  End  time  scEles  of  the 
secondEry  flows  under  considerEtion,  the  pETEllel  instEbility  mechEnism  is  of  little 
consequence  Et  our  scEle  (Brown  1980).  The  third  mechEnism,  thermEl  instEbility,  derives 
its  forcing  fi'om  the  Eir-seE  temperEture  difference  snd  CEn  be  likened  to  clESsic  Rnyleigh- 
BenErd  convection.  For  consistency  with  Shirer  et  al.  (1996),  we  refer  to  the  rolls  End 
cells  ES  boundEry  kyer  convective  eddies  (BLCE).  The  verticEl  extent  of  the  BLCE  is 
given  by  the  boundEry  kyer  depth.  End  their  primEry  effect  is  the  stEbilizEtion  of  the 
Etmosphere  through  the  trEnsport  of  heEt  End  momentum. 
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Numerous  projects  have  been  undertaken  to  model  BLCE  using  nonlinear  dynamical 
systems  (e.g.,  Shirer  1986;  Haack  and  Shirer  1992;  Laufersweiler  and  Shirer  1995).  The 
primary  goal  of  these  studies  was  to  determine  the  temporal  and  spatial  responses  of 
BLCE  to  different  forcing  rates,  both  mechanical  and  thermal.  These  models  successfully 
capture  rolls  and  provide  the  fluxes  of  heat  and  horizontal  momentum  as  well.  However, 
due  to  the  limitations  imposed  by  the  rigid,  stress-free  and  perfectly  conducting  upper  and 
lower  boundary  conditions,  we  are  unable  to  use  any  of  these  models  to  deduce  the  role 
BLCE  play  in  perturbing  the  sea  surface,  as  is  required  to  apply  the  remote  sensing 
capabilities  of  SAR  to  the  MABL. 

To  overcome  this  obstacle,  Zuccarello  (1994)  developed,  and  Lambert  (1995)  later 
improved,  a  nonlinear,  three-dimensional,  Boussinesq  model  that  provides  for  subgrid- 
scale  momentum  and  heat  fluxes  at  the  lower  boundary.  Two  constant  forcing  parameters 
are  derived  from  Monin-Obukhov  similarity  theory  to  capture  the  flux  contributions  in  the 
surface  layer;  these  yield  thermodynamic  forcing  and  momentum  dissipation.  The  new 
boundary  conditions  are  more  realistic  and  can  be  used  to  deduce  wind-imposed  stress  on 
the  sea  surface.  Following  the  success  of  Laufersweiler  and  Shirer  (1995),  Lambert 
introduced  (after  his  1995  work)  an  inversion  layer  into  the  model  to  cap  moderately 
forced  secondary  circulations.  However,  we  still  specify  the  upper  boimdary  to  be  rigid, 
stress-free  and  perfectly  conducting,  as  was  done  in  previous  MABL  models. 

The  initial  results  obtained  by  Zuccarello  (1994)  in  his  case  study  were  encouraging, 
for  they  indicated  that  the  model  adequately  captures  the  spatial  characteristics  of  the 
MABL  secondary  circulations.  Additionally,  the  component  of  the  mean  wind 
perpendicular  to  the  roll  axis  (i.e.,  the  cross-roll  mean  wind)  and  the  period  of  the  system 
reflect  those  that  are  typically  observed  in  the  atmosphere.  However,  Zuccarello  did 
encounter  several  problems  during  his  case  study.  One  was  the  presence  of  an  unknown 
energy  source  that  prevented  model  equilibration  during  integration.  The  second  was  a 
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vertical  velocity  maximum  that  appeared  at  the  lower  boundary,  instead  of  within  the 
boundary  layer  at  the  typically  reported  height  of  0.3  times  the  domain  depth  (Lenschow 
et  al.  1980).  Lambert  (1995)  conducted  an  energetics  analysis  and  identified  an 
extraneous  kinetic  energy  source  at  the  lower  boundary  of  the  model.  By  imposing  two 
additional  boundary  conditions  on  the  flow,  he  eliminated  the  unwanted  energy  source  as 
well  as  the  vertical  velocity  anomaly,  enabling  model  integration  and  quasi-equilibration. 
Lambert  also  developed  a  subroutine  to  construct  the  model  output  and  aid  in  the 
visualization  of  the  model  solution.  The  apex  of  his  work  was  the  creation  of  an  x-y 
planview  that  depicts  reasonable  wind-imposed  stress  patterns  on  the  sea  surface. 

The  purpose  of  this  study  is  to  continue  Zuccarello  (1994)  and  Lambert’s  (1995)  work 
on  the  impact  BLCE  have  on  sea  surface  stress,  with  an  emphasis  on  visualization  of  the 
model  solutions.  In  Chapter  2  we  summarize  background  theory  and  model  development, 
then  re\iew  the  system  of  partial  differential  equations  and  the  spectral  conversion.  We 
conclude  the  chapter  by  describing  the  details  of  the  computer  model.  In  Chapter  3  we 
note  improvements  made  to  the  model  since  Lambert’s  (1995)  work,  and  in  Chapter  4,  we 
present  case  studies  from  the  Hi-Res  2  field  experiment  (Lambert  1995;  Sikora  et  al 
1995).  We  conclude  by  proposing  areas  for  future  research. 
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Chapter  2 

MODEL  OVERVIEW 

> 

Boundary  layer  convective  eddies  (BLCE)  have  been  the  focus  of  many  studies  in 
recent  years.  A  number  of  these  employ  nonlinear  dynamical  systems  to  model  kilometer- 
scale  circulations  (e.g.,  Shirer  1986;  Haack  and  Shirer  1992;  Laufersweiler  and  Shirer 
1995).  The  motivation  behind  these  studies  has  been  to  explain  the  temporal  and  spatial 
responses  of  BLCE  to  thermal  and  mechanical  forcing  mechanisms.  Such  models 
successfully  capture  the  secondary  circulations  and  provide  the  associated  vertical  flux 
profiles  of  heat  and  horizontal  momentum-but  within  the  boundary  layer,  not  across  the 
boundaries  themselves.  For  simplicity,  the  upper  and  lower  boundaries  used  in  these 
models  are  rigid,  stress-free  and  perfectly  conducting.  Hence,  fluxes  at  or  near  the 
boundaries  can  not  be  adequately  represented  and  the  effects  of  BLCE  on  sea  surface 
stress  can  not  be  deduced,  as  is  required  for  the  application  of  SAR  technology  to  the 
study  of  the  MABL. 

Zuccarello  (1994)  develops  a  spectral  model  that  incorporates  nonzero  subgrid-scale 
momentum  and  heat  fluxes  at  the  lower  boundary,  allowing  the  study  of  BLCE-induced 
sea-surface  stress.  In  specifying  the  boundary  conditions,  he  invokes  similarity  theory 
(Stull  1988).  For  simplicity  and  consistency  with  earlier  models,  Zuccarello  keeps  the 
upper  vertical  boundary  rigid,  stress-firee  and  perfectly  conducting.  Lambert  (1995) 
expands  upon  Zuccarello’s  work  in  several  areas:  he  (a)  conducted  an  energetics  analysis 
to  identify  an  extraneous  energy  source  and  vertical  velocity  anomaly  at  the  lower 
boundary;  (b)  reformulated  the  lower  boundary  conditions  to  eliminate  the  problems  in  (a); 
and  (c)  developed  an  output  subroutine  to  aid  the  analysis  and  interpretation  of  the 
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model  solution,  with  an  emphasis  on  the  lower  boundary  stress  planview  for  comparison 
with  SAR.  After  finishing  his  thesis,  Lambert  incorporated  an  inversion  subroutine  into 
the  computer  model. 

> 

Here,  we  extend  the  work  done  by  Zuccarello  and  Lambert  on  their  BLCE  model, 
focusing  primarily  on  visualization  of  the  model  solutions. 

2. 1 .  Summary  of  Model  Development  and  Background  Theory 

Zuccarello  (1994)  provides  a  detailed  account  of  model  development  and  background 
theory.  Here,  we  summarize  Zuccarello’s  work,  as  did  Lambert  (1995),  and  encourage 
interested  readers  to  review  his  unabridged  presentation.  We  have  modified  Zuccarello 
and  Lambert’s  temperature  base  state  to  allow  for  the  inclusion  of  an  inversion  and  an 
unstable  surface  layer  of  the  tj^ie  used  by  Laufersweiler  and  Shirer  (1995). 

We  employ  the  shallow  Boussinesq  system  of  equations  to  capture  the  convective 
nature  of  boundary  layer  secondary  circulations.  These  circulations,  manifested  as  two- 
dimensional  rolls  or  three-dimensional  cells,  generally  transport  heat  upward  and 
horizontal  momentum  downward,  while  tending  to  alter  the  vertical  shear  in  the  boundary 
layer. 

The  basic,  or  conductive  state,  is  considered  to  be  time-independent,  isopycnic, 
hydrostatic,  incompressible  and  horizontally  moving  (Shirer  1986).  It  is  represented  by 
the  standard  meteorological  variables  with  subscript  0  in  the  following  equations; 


v.(z)  =  V(z)  =  £/(r)i+^'(z)j 


Pq  ~  Poo 


(2.1) 

(2,2) 
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^o(^)  =  ^oo  +  ^.w  (2-3) 

aW  =  /^oo-Poo^^>  (2-4) 

<  where  the  00  subscript  denotes  atmospheric  values  at  the  bottom  of  the  domain  and 

is  the  height-dependent  portion  of  the  temperature  profile.  We  have  rewritten  the  thermal 
base  state  using  potential  temperature  9q{z)  ,  which  is  composed  (in  the  model)  of  three 
different  linear  profiles  in  the  surface  layer,  the  neutral  boundary  layer  and  the  inversion 
layer  (Section  3.1).  Moreover,  the  mean  MABL  density  is  essentially  equal  to  the  density 
at  the  sea  surface  (2.2). 

We  represent  rolls  and  cells  as  perturbations  superimposed  on  this  basic  state. 
Variables  unique  to  the  secondary  circulation  perturbations  are  denoted  by  primes  in  the 


following  equations; 

u'{x,y,z,t)  =  t{x,y,z,t)-U{z)  (2.5) 

v'{x,y,z,t)  =  v{x,y,z,t)-V(z)  (2.6) 

w'{x,y,z,t)  =  w{x,y,z,t)  (2.7) 

ff[x,y,z,t)  =  f:{x,y,z,t)-p^  (2.8) 

e\x,y,z,t)  =  e(x,y,z,t)  -  e^{z)  (2.9) 

p\x,y,z,t)  =  p{x,y,z,t)  -  pj^z).  (2. 10) 


We  have  chosen  the  horizontal  domain  to  be  infinite  and  cyclically  continuous,  with 
gravest  modes  Lx  and  Ly,  where  Lx  and  Ly  are  the  x  andj  domain  lengths  necessary  to 
encompass  one  entire  wavelength  L  of  the  circulation  =  Lx^  +  Ly^).  The  vertical 

domain  extends  fi'om  the  lower  boundary,  hLB=\0m  (anemometer  height),  to  the  cloud 
top  height  (zr).  Figure  2.1,  which  shows  one  half  of  the  horizontal  domain,  clearly  depicts 
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the  secondary  circulation  as  bounded  by  the  base  of  the  inversion  z,  and  hiB-  Since 
temperature  decreases  and  wind  speed  increases  with  height,  we  can  visualize  in  the 
air/sea  interface  layer  how  the  net  vertical  heat  flux  at  z  =  hiB  should  be  positive  and  the 
net  vertical  transport  of  horizontal  momentum  should  be  negative. 

We  rewrite  the  entire  domain  in  dimensionless  form  as  (0,0,0)  <Jx,y,z')  <.(1,1,1). 
The  model  is  designed  to  capture  the  circulation  above  hiB  and  employs  surface  layer 
similarity  theory  (Stull  1988)  to  parameterize  the  transport  of  heat  and  momentum  below 
via  the  boundary  conditions  described  in  the  next  section. 


MODEL  BLCE  CIRCULATION 


Fig.  2. 1.  A  schematic  cross  section  of  a  BLCE  in  the  MABL,  in  which  the  air/ocean  interface  layer  is 

greatly  exaggerated.  The  air/ocean  interface  layer  potential  temperature  decreases  and  the  wind 
speed  increases  with  height.  The  model  domain  extends  from  z  =  hiB  (z  =  0)  to  z  =  zr  (z  =  1) 
(Adapted  from  Shirer  et  al.  1996). 
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2.2.  The  System  of  Partial  Differential  Equations 

In  order  to  make  our  Boussinesq  system  of  equations  dimensionless,  Zuccarello 
(1994)  offers  the  following  definitions; 


V'  - 

(2.11) 

,  W*  K 

^  ~z  -A 

(2.12) 

x  =  x*L^ 

(2.13) 

• 

1! 

(2.14) 

z  —  z  {zj 

(2.15) 

K 

(2.16) 

Q,  _ 

g(zj, 

(2.17) 

(zj.  “* 

(2.18) 

VW  =  |V(r,)|vy) 

(2.19) 

/= 

(2.20) 

d  ,  d  .  d  . 

(2.21) 

where  the  constant  thermal  conductivity  is  represented  by  k,  the  constant  eddy  viscosity 
by  V  and  the  roll  aspect  ratios  by  and  Oy,  defined  below.  We  then  eliminate  pressure 
from  the  momentum  equation  by  taking  the  curl  to  form  a  vorticity  equation.  Together 
wth  the  thermodynamic  and  continuity  equations,  the  resulting  dimensionless  system  of 
equations  is 
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^  =  V'6>*  -  /?eV’  •  V^*  -  -  V*  • 

dt  a. 

=  ^V'(V  X  V*)  -/*P[V  X  (k  X  V*)]  +  P(V  X  ^k)  -  i?e[v  x  (v*  •  Vv*)] 


(2.22) 


-^V  X  |v*  •  Vv*|j  -  Re 


Vx  w  — 
V«v*  =  0. 


(2.23) 


(2.24) 


In  (2.22),  Lambert  (1995)  chose -(<^o7^*)  =  -^^  =  “Ito  represent  the  background 
neutral  boundary  layer.  We  can  therefore  extract  five  dimensionless  parameters  from  our 
system  of  equations-the  eddy  Prandtl  number,  P;  the  roll  aspect  ratios,  and  ;  the 

Reynolds  number.  Re;  and  the  Rayleigh  number,  Ra.  Their  definitions  are 


P  =  vIk 

(2.25) 

II 

(2.26) 

_  {^T  ^Lb) 

L 

{2.21) 

Re  =  |v(z2-)|(zr 

(2.28) 

~^Lb) 

(2.29) 

vkT^ 

where  is  the  potential  temperature  difference  across  the  boundary  layer.  The 
Reynolds  number  characterizes  the  dynamic  forcing  by  the  mean  wind  shear  in  the  basic 
state,  while  the  Rayleigh  number  represents  thermodynamic  forcing  throughout  the  (nearly 
neutral)  boundary  layer. 
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> 


For  the  dimensionless  Boussinesq  system  of  differential  equations,  we  choose  the 
horizontal  domain  to  be  cyclic  and  specify  the  upper  vertical  boundary  to  be  rigid,  stress- 
free  and  perfectly  conducting.  The  lower  boundary,  however,  is  an  important  source  of 
thermal  energy  for  secondary  circulation  development  and  a  sink  for  momentum  (Fig.  2.1). 
In  specifying  lower  boundary  conditions,  Zuccarello  (1994)  defines  two  parameters  to 
control  the  forcing:  the  Schramm  momentum  constant,  Sm,  and  the  Schramm  temperature 
constant  st,  both  of  which  are  positive,  consistent  with  Figure  2. 1 .  The  upper  and  lower 
boundary  conditions  are 


o 

II 

(2.30) 

o 

II 

• 

II 

* 

II 

(2.31) 

«'(o)+^^=o 

(2.32) 

(2.33) 

(2.34) 

^•(0)  _  <?V(0) 
az*  ”  dz^ 

(2.35) 

Note  that  (2.35),  which  follows  from  (2.32)  and  (2.34)  via  the  continuity  equation,  is  not 
ideal  in  that  it  involves  derivatives  of  the  vertical  velocity  w  at  the  lower  boundary, 
instead  of  w*  itself.  In  early  model  runs,  this  anomaly  manifested  itself  by  producing 
vertical  velocity  maxima  at  the  lower  boundary,  instead  of  the  more  physically  reasonable 
w*  =  0.  To  overcome  this  obstacle,  Lambert  (1995)  introduces  two  new  boundary 
conditions: 


w*(0)  =  0 
A*(0)  =  0, 
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(2.36) 

(2.37) 

where  A*  is  the  vector  streamfiinction  defined  in  the  following  section.  In  the  atmosphere, 
the  BLCE  vertical  velocity  perturbation  at  /ilb  =  10  w  (z*  =  0)  is  not  identically  zero,  but 
to  an  excellent  approximation  we  may  regard  it  as  such. 

To  attain  boundary  conditions  (2.36)  and  (2.37),  Lambert  (1995)  redefines  the  vertical 
streamfiinction  basis  functions  instead  of  rewriting  the  entire  model.  For  a  complete 
account  of  this  procedure,  we  refer  the  reader  to  Appendix  B  of  his  work. 

2.3.  From  Differential  Equations  to  Computer  Code:  The  Spectral  Conversion 

To  create  a  spectral  model  from  our  system  of  partial  differential  equations  (PDEs), 
Zuccarello  (1994)  employs  the  Galerkin  technique  or  spectral  method  (Higgins  1987). 
When  applied  to  the  vorticity  (2.23)  and  thermodynamic  (2.22)  equations,  this  method 
yields  a  nonlinear  dynamical  system  (NDS).  Originally,  the  NDS  consisted  of  60 
equations--20  for  temperature  and  40  for  streamfiinction-but  was  reduced  to  50  when 
Lambert  (1995)  introduced  new  lower  boundary  conditions  (2.36)  and  (2.37).  The  new 
NDS  still  has  20  equations  for  temperature  but  now  has  only  30  for  streamfiinction. 

Since  the  velocity  perturbation  vector  v*  is  nondivergent,  it  can  be  represented  by  the 
curl  of  a  vector  streamfiinction  A*,  whose  components  are  y/  and  tj  in  the  x  -  and^  - 
directions,  respectively: 


v*  =  VxA*  =Vx(</*i)+Vx(;7*j). 


(2.38) 
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Truncated  Fourier  expansions  for  9',  y/  and  rf  are  composed  of  time-dependent 
amplitudes,  spatially-dependent  trigonometric  functions  with  horizontal  wavenumbers  / 
and  m,  and  vertical  basis  functions,  F//)  and  h^z*): 


4  4 


(r(x\y\2/)=j;L^,(t-)trig,(x\yy,(z) 

1=0  y=i 

(2.39) 

W(x\y,z/)  =  ttw„(>-)>ngA^'yW) 

p=0  9=l 

(2.40) 

=  ±t  ’?„{<yig,(x.y)h,{z-) , 

p=0  q=\ 

(2.41) 

where 

trigo(x\y)  =  1 

(2.42) 

trigi[x*,y*^  =  sin2;^x*  smlTWty* 

(2.43) 

trigj^x*,y*)  =  sin2;^x*cos2;z7wj* 

(2.44) 

trigj^x*  ,y*^  =  co&lrdx*  smlTony* 

(2.45) 

trigj^x*,y*^  -  QoslTdx*  coslnmy* . 

(2.46) 

Here  the  orthogonal  ^‘vertical  basis  functions  are 

F.{z*\  =  cosh^iZ* - ^sinho.z*  for  5j.  <  1 

(2.47) 

F,{z*\  =  cos (o.z* - —  sin 6),z*  for  Jj.  >  1 , 

(2.48) 

where 

SjCO^  =  tanhOj  for  5^.  <  1 

(2.49) 

=tanfUi  for  5j.>l, 

(2.50) 

and 
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1  . 


F, (z*)  =  COSO ,z* - sin o  z*  for  j  =  2,3,4 , 

J  \  /  J  c*  /i-fc  ^ 


(2.51) 


where 


Sj.(Oj  =  tan  cOj  for  j  =  2,3,4 . 


(2.52) 


The  nonorthogonal  t//  and  rf  vertical  basis  functions  are  expressed  as  differences: 


/»,(z*)  =  (cos  C7,z*  -  s„mg  sin  ar^z*)  -  (cos  sin  nr^+,z*)  for  ^  =  1, 2, 3 ,  (2.53) 

where 

^  ~  1>2,3,4.  (2.54) 

To  produce  our  NDS  (2.39)  -  (2.41)  using  the  standard  Galerkin  procedure  (Higgins 
1987),  we  insert  the  0\  y/  and  t]  expansions  into  the  dimensionless  PDEs  multiplied  by 
the  corresponding  basis  functions,  and  integrate  over  the  domain  (Shirer  1987). 
Appropriate  trigonometric  identities  are  applied  to  simplify  the  resulting  terms  in  the 
equations. 

In  the  NDS,  each  PDE  term  is  represented  as  the  product  of  a  50  X  50  matrix  and  a 
50-element  amplitude  coefficient  vector,  allowing  us  to  write  the  model  equations  in  the 
form: 

Ay  =  By  +  Ciy)y,  (2.55) 

where  .<4  is  a  50  X  50  invertible  matrix  of  inner  products;  5  is  a  50  X  50  constant  linear 
matrix,  determined  by  P,  Re,  Ra,  a, ,  and  the  Fourier  coefficients  of  the  background 

temperature  and  wind;  and  C  is  a  50  X  50  linear-in-j^  matrix  from  the  advective  terms  of 
the  two  systems  of  PDEs.  The  50-element  column  vector;;  contains  20  time-dependent 
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spectral  amplitudes  for  temperature,  TjJ ,  and  30  time-dependent  spectral  amplitudes  for  the 
streamfimction  components,  and  7^, .  The  50-element  column  vector  y  represents 

the  temporal  derivatives  of  the  spectral  amplitudes  in_y.  Since  A  is  invertible,  we  may 
rewrite  (2.55)  in  the  form 

y^A-'By+A-^C{y)y.  (2.56) 

This  is  the  form  used  in  the  numerical  integration. 

2.4.  The  Computer  Model 

The  physical  and  mathematical  model  described  in  the  preceding  sections  was  encoded 
by  Zuccarello  (1994)  using  the  FORTRAN  programming  language;  it  was  originally 
compiled  with  a  UNIX-based  WATFOR77  compiler.  To  add  portability,  flexibility  and 
ease  of  operation,  we  developed  a  PC  version  of  the  code  using  Microsoft’s  Fortran 
Power  Station  4.0  software.  Code  editing,  compilation  and  execution  are  now  done 
easily. 

The  sequence  of  steps  necessary  to  carry  out  a  SAR  analysis  using  the  model  are 
summarized  by  the  flowchart  in  Figure  2.2  and  are  referenced  periodically  throughout  the 
text  that  follows.  As  with  any  study,  we  first  select  a  case  (Figure  2.2,  Step  A)  and 
specify  atmospheric  input  parameters  like  those  shown  in  Table  4.1. 

Using  a  subroutine  discussed  in  Appendix  A  of  Zuccarello  (1994),  we  input  the  wind 
speed  and  mixing  ratio  at  10  /w  as  well  as  the  air-sea  temperature  difference,  and  employ 
Monin-Obukhov  surface  layer  similarity  theory  to  determine  values  for  the  Schramm 
momentum  and  temperature  constants,  Sm  and  st  (Figure  2.2,  Step  B).  Based  upon  these 


A  SAR  ANALYSIS 
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Fig.  2.2.  Summary  of  required  steps  in  a  SAR  analysis. 
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input  values,  the  code  calculates  values  for  ©j  and  tUc,  where  j,c  =  1,2, 3, 4,  and  then 
iteratively  solves  (2.49),  (2.50),  (2.52)  and  (2.54)  for  the  eigenvalues.  Lambert  (1995) 
empirically  derives  a  single  set  of  initial  oj  and  etc  values  for  this  iterative  process  that 
works  well  for  most  cases.  However,  for  certain  values  of  s„  and  sr,  we  encounter 
problems,  suggesting  that  ©j  and  ©care  jumping  to  other  quadrants  of  their  respective 
tangent  and  cotangent  fiinctions;  tWs  phenomenon  causes  an  undesirable  convergence 
to  previously  determined  omega  values.  To  correct  this  problem,  we  choose  seed  values 
that  are  multiples  of  k  for  qj,  and  multiples  of  7c/2  for  TOc;  for  the  case  sr  <  1,  we  set  = 
(srXK  Then  we  apply  the  bisection  method  to  approximate  each  ©j,  ©c,  thus  obtaining 
initial  values  for  the  Newton’s  method  iteration  scheme.  Extensive  testing  has 
demonstrated  the  reliability  of  this  automated  eigenvalue  routine. 

We  next  construct  a  background  wind  profile  using  the  four  fi'ee  parameters;  (/(zr), 
V(zr),  du  and  dy  (Figure  2.2,  Step  C).  The  former  two  are  chosen  by  the  user  to  control 
the  overall  speed  and  turning  shear  in  the  profile,  while  the  latter  two  are  the  free 
constants  in  the  wind  profile  equations: 

U(z)  =  a^jZ+hjjZ^'^  +Cjjz''^  +d^j  (2.57) 

V(z)  =  ayZ-¥byz'^'^  +<Vz''^  +dy .  (2.58) 

We  apply  (2.15)  and  (2.19)  to  nondimensionalize  these  profile  equations  and  then  require 
the  mean  background  wind  profiles  to  satisfy  the  same  boundary  conditions  as  the 
perturbation  wind  velocities  do.  These  boundary  conditions  are 


./  .  \  ^*{z*  =  0) 


(2.59) 
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3^*(z*  =  O) 

vy  =  0)-s„  . i  =  0.  (2.60) 

Additionally,  we  must  specify  C/(z*  =  0)  and  V(z*  =  0),  which  are  the  same  boundary  values 
used  to  calculate  the  Schramm  momentum  constant  Sm.  As  a  result  of  our 
nondimensionalization  definition  (2.19),  we  note  that  the  following  condition  is  also  met; 

|v*(z*  =  l)|  =  l.  (2.61) 

Because  (2.61)  is  not  a  new  boundary  condition,  we  are  left  with  she  equations  and  eight 
unknown  variables,  two  of  which  are  the  fi'ee  wind  parameters  du  and  dy.  Lambert  (1995) 
developed  an  EXCEL  spreadsheet  to  enable  users  to  easily  find  the  values  for  du  and  dy 
necessary  to  give  physically  realistic  background  wind  profiles.  Figure  2.3  shows  an 
example  of  how  the  wind  profile  speed  changes  by  altering  the  values  oidu  and  dy.  For 
this  scenario,  we  experiment  until  we  arrive  at  the  best  possible  wind  profile  for  the  given 
atmospheric  conditions-curve  B  in  Figure  2.3.  This  plausible  boundary  layer  background 
wind  profile  corresponds  to  c?£/  =  4. 1  and  dy=\.5  when  |V(hLB)|  =  4  m/s  and  |V(zt)1  =  8 
m/s.  Wind  direction  profiles,  though  not  shown,  are  determined  concurrently  with  the 
speed  profiles.  Before  moving  fi-om  Step  C  to  Step  D  in  Figure  2.2,  we  must  specify  a 
plausible  inversion  strength  (®C/m),  preferably  from  an  atmospheric  sounding. 

A  linear  stability  analysis  (Figure  2.2,  Step  D)  is  then  performed  to  determine  what 
domain  shape,  specified  by  Lx  and  Ly,  yields  the  largest  growth  rates  for  a  given  domain 
depth  zp— these  are  the  modes  we  would  expect  to  observe  in  the  atmosphere 
(Laufersweiler  1987).  With  the  assistance  of  a  contouring  software  package,  we  can 
readily  identify  the  fastest  growing  mode  and  its  corresponding  and  Oj,  values,  the 

preferred  aspect  ratios.  Before  integrating  our  system  of  equations,  we  determine  a  set  of 
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BACKGROUND  WINDSPEED 


du  =  3.1;  dv  =  0.5  (A) 
du  =  4.1;  dv  =  1.5  (B) 
du  =  5.1;  dv  =  2.5  (C) 


0.5  0.75  1.0  1.25 


Fig.  2.3.  Dimensionless  wind  speed  profiles  as  determined  by  the  mean  background  wind  parameters,  du 
and  dv-  Note  how  the  wind  speed  profile  changes  from  A  to  C  as  du  and  dy  are  incremented  by  1.0.  We 
find  du  and  dy  empirically  to  attain  more  physically  plausible  wind  profiles.  Here,  we  choose  profile  B. 
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initial  conditions  for  the  50-element  vector^  using  the  eigenvectors  associated  with  the 
eigenvalues  ©j  and  ©c  (Figure  2.2,  Step  F).  If  we  desire  a  continuation  run,  then  we  use^'- 
vector  values  from  the  end  of  the  previous  integration. 

Matrix^  and  matrix 5  (2.55),  which  represent  the  temporal-derivative  terms  and  all 
other  linear  terms,  respectively,  are  determined  prior  to  integrating  because  they  do  not 
depend  on  the  spectral  coefficients  In  contrast,  the  nonlinear  term  matrix  C  is  a  linear 
function  ofy  and  must  be  recomputed  at  each  integration  time  step,  at  least  in  theory. 
With  efficiency  in  mind,  Lambert  (1995)  re-expresses  the  linear  function  C(y  )  using  a 
cubic  tensor  of  coefficients  (CCUBE,  Section  3.2)  so  that  these  need  not  be  recomputed 
at  each  time  step.  Of  course,  the  same  shortcut  applies  to  the  coefficient  term  A'^C(y  )  in 
(2.56).  We  use  the  IMSL  subroutine  DLINRG  to  convert  matrix  A  into  A'^,  employ 
DMRRRR  to  determine  A'^B,  and  calculate  A'^C  directly. 

Earlier  versions  of  the  model  utilized  the  IMSL  subroutine  DIVPAG  to  perform  the 
integration  (Figure  2.2,  Step  G),  with  acceptable  results  for  short  (i.e.,  less  than  32,000 
time  step)  runs.  For  longer  runs,  however,  we  encountered  eigenvalue  convergence 
problems  and  numerical  noise,  suggesting  that  another  integration  scheme  is  necessary 
(Section  3.2).  Since  our  system  of  equations  is  stif^  we  adopted  the  semi-implicit 
integration  scheme  STIFBS  from  Press  et  al.  (1992).  Though  not  as  fast  as  DIVPAG, 
STIFBS  has  greatly  abated  the  noise  and  convergence  problems,  yielding  smooth, 
consistent  integrations. 

We  allow  the  model  to  integrate  until  the  dimensionless  energy  {y^ Ay  )  quasi- 
equilibrates.  After  the  last  integration  time  step,  the  model  passes  the  time-dependent 
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amplitude  coeflBcients  to  a  computation/output  (CO)  subroutine  developed  by  Lambert 
(1995),  where  the  solution  is  constructed  using  (2.38)-(2.41).  We  have  enhanced 
Lambert’s  CO  subroutine  by  offering  additional,  more  automated  options. 

In  Chapter  4,  we  offer  examples  of  model  solutions  generated  by  the  CO  subroutine 
for  the  Hi-Res  2  field  experiment.  Before  doing  so,  however,  we  proceed  to  the  next 
chapter  where  we  expound  the  improvements  made  to  the  model  since  Lambert  (1995). 
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Chapter  3 

MODEL  IMPROVEMENTS 


At  the  onset  of  this  study,  we  had  hoped  to  immediately  resume  work  on  the 
interpretation  of  a  SAR  image  taken  during  Hi-Res  2,  a  1993  MABL  field  experiment 
conducted  off  the  coast  of  North  Carolina  (Lambert  1995;  Sikora  et  al.  1995).  However, 
early  model  runs  using  Hi-Res  2  observations  as  input  parameters  were  quite 
discouraging.  In  addition  to  very  high  frequency  noise  problems  in  the  energy, 
streamfunction  and  temperature  perturbations,  we  also  observed  startling  dimensionless 
energy  growth  rates  that  terminated  model  runs  prematurely.  We  experimented  with 
altered  initial  conditions  and  subroutines  to  limit  numerical  roundoff  errors  with  little  or  no 
effect,  convincing  us  that  further  model  improvements  and  debugging  were  necessary.  We 
also  were  compelled  to  address  previously  unexplained  numerical  issues  related  to  our 
NDS.  Therefore,  before  pursuing  the  Hi-Res  2  case  study  further,  we  had  to  address 
model  limitations  and  make  corrections  and  modifications  as  necessary. 

Preliminary  troubleshooting  revealed  that  the  immense  growth  rates  were  due  to  the 
terms  represented  by  the  linear  buoyancy  matrix,  and  that  the  noise  problem  was  generated 
by  several  sources,  two  of  which  were  1)  the  integration  of  stiff  equations  with  a  non-stiff 
numerical  integration  routine  (Press  et  al.  1992),  and  2)  the  nonlinear  cubic  tensor  arising 
fi-om  the  advection  of  vorticity  (Section  2.4).  Since  this  growth  rate  problem  generally 
caused  the  computer  code  to  crash,  we  decided  to  focus  on  the  linear  portion  of  the  model 
first. 
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3.1.  Addressing  Monotonic  Energy  Growth 

While  constructing  the  model,  Zuccarello  (1994)  built  all  the  linear  matrices 
analytically  using  the  DERIVE  software  package.  After  completing  his  1995  case  study, 
Lambert  did  likewise  when  he  constructed  the  inversion  subroutine,  the  terms  of  which  he 
merely  added  to  those  representing  the  neutral  boundary  layer  (i.e.,  Ra  =  -1),  as  shown  in 
the  equation 

^  =  -RaA.  (3.1) 

a  oz 

In  order  to  perform  a  completely  independent  cross-check  of  the  linear  buoyancy  matrix, 
we  began  with  the  spectral  representations  for  9  *,  y/  and  t]  (2.39)-(2.41),  encoded  our 
own  fimctions  and  integrated  the  results  numerically.  Since  the  trigonometric  functions 
(2.42)-(2.46)  are  independent  ofz*,  and  the  basis  functions  (2.47)-(2.52)  are  independent 
of  x*  andy,  we  were  able  to  compute  the  x*,  y,  z  triple  integral  as  a  product  of  a  double 
integral  in  x  ,  y  and  a  single  integral  in  z  .  We  used  the  IMSL  subroutine  DTWODQ  for 
the  double  integral  and  DQD  AG  for  the  single  one. 

The  analytic  and  numeric  approaches  to  calculating  the  inversion  were  nearly  identical, 
indicating  that  the  inversion  subroutine  had  been  coded  correctly.  However,  a  closer 
examination  of  the  analytically  derived  neutral  boundary  layer  revealed  unexplained  errors. 
In  the  interest  of  time,  we  chose  to  incorporate  the  numerical  subroutine  into  the  model 
rather  than  to  reformulate  the  analytic  equations.  This  decision  had  the  added  benefit  of 
making  future  modifications  of  the  temperature  profile  more  easy  to  implement. 

After  fiilly  incorporating  the  numerically  determined  temperature  profile,  however,  we 
found  that  the  introduction  of  an  inversion  still  tended  to  increase  the  overall  energy  of  the 
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system,  even  with  very  weak  inversions.  The  culprit  appeared  to  be  the  dimensionless 
inversion  lapse  rate  term 

^  =  (3-2) 

d:  & 

whose  immense  magnitude  (of  order  10*)  dominated  all  other  terms  in  the  NDS,  especially 
those  of  the  background  temperature  modification  9^^  (2.39).  Adding  an  inversion  to  a 

thermally  neutral  boundary  layer  (of  order  unity)  tended  to  throw  the  model  solution  out 
of  balance,  yielding  colossal  growth  rates  in  the  dimensionless  energy. 

Hoping  to  counteract  this  aphysical  growth  and  yield  more  balanced  solutions,  we 
introduced  an  unstable  surface  layer  that  spreads  the  air-sea  temperature  difference  over  a 
user-specified  surface  layer  depth  (Laufersweiler  and  Shirer  1995);  this  yields  a  lapse  rate 
term  of  similar  magnitude  to  (3.2).  The  new  potential  temperature  profile  is  similar  to  the 
one  depicted  in  Figure  3.1,  although  the  depth  of  the  unstable  surface  layer  has  been 
exaggerated  somewhat.  The  three  lapse  rates  are,  by  definition,  negative  in  the  surface 
layer,  slightly  positive  in  the  boundary  layer  and  positive  in  the  inversion  layer. 
Quantitatively,  the  dimensional  potential  temperature  lapse  rates  can  be  represented  by  the 
following  equations: 


^SL 

(hLB  <  Z  <  Zsl) 

(3.3) 

dz 

^BL 

< 

1 

{ZSL  <  Z  <  Zi} 

(3.4) 

& 

dOj^ 

<1 

1 

{Zi  <  Z  <  Zt}, 

(3.5) 

& 

HEIGHT 
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THERMAL  BASE  STATE 


POTENTIAL  TEMPERATURE 


Figure  3.1.  A  schematic  of  the  background  potential  temperature  profile  used  in  the  model. 
(Not  drawn  to  scale.) 
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where  zsl  is  the  top  of  the  unstable  surface  layer  and  is  the  change  in  potential 
temperature  over  the  specified  layer,  n.  All  other  variables  are  as  defined  in  Chapter  2. 
This  addition  to  the  model  did  increase  the  magnitude  of  the  other  buoyancy  terms  as 
desired  but  did  not  completely  balance  the  model  solutions  as  we  had  hoped—there  was 
still  evidence  of  monotonic  energy  growth  in  several  instances.  We  therefore  added  a 
nonlinear  numerical  dissipation  term  to  bound  the  energy  growth,  without  altering  the 
linear  growth  rates  or  the  preferred  aspect  ratios.  Equation  (2.56)  now  takes  the  form: 

y^A-^By^A-^Ciy)y-d\y\y,  (3.6) 

where  d  is  an  empirically  determined  dissipation  coefficient;  we  use  d=0.\  for  our  study. 
Model  runs  with  an  inversion  and  surface  layer  now  provide  more  stable,  physically 
plausible  results. 

3.2.  Addressing  Numerical  Noise 

Isolating  the  source  of  numerical  noise  in  the  nonlinear  terms  proved  to  be  less 
straightforward  because  the  problem  did  not  occur  all  the  time.  Model  output  showing 
the  emergence  of  noise  was  extremely  dependent  on  initial  conditions,  leading  us  to 
consider  the  possibility  that  the  solution  was  chaotic  (Lorenz  1993).  As  it  turned  out, 
such  was  not  the  case. 

In  the  computer  model  (Section  2.4),  there  is  a  50  X  50  X  50  matrix,  CCUBE,  that  is 
used  to  calculate  the  C(y)y  term  in  (2.56).  This  cubic  matrix  holds  the  constant 
coefficients  for  the  nonlinear  temperature  and  streamfunction  advective  terms.  Numerous 
diagnostic  tests  revealed  that  several  CCUBE  assignment  statements  contained  indexing 
errors  that  manifested  themselves  as  noise  in  the  model  solution.  After  Drs.  Hampton 
Shirer  and  Robert  Wells  tediously  re-indexed  the  CCUBE  matrix,  we  noted  vast 
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improvements  in  the  model  solution;  the  nonlinear  terms  became  energy  preserving  as 
they  should  and  numerical  noise  was  minimized,  occurring  less  frequently.  However, 
these  changes  did  not  totally  eradicate  numerical  noise  from  the  model  output  as  we 
would  have  liked,  for  during  longer  integrations  when  surface  winds  were  light,  noise 
continued  to  dominate  the  solution.  The  installation  of  a  time  filter  reduced  the  amplitude 
of  the  fluctuations  somewhat  but  did  not  remove  the  higher-frequency  oscillations. 

As  noted  in  Chapter  2,  we  encountered  eigenvalue  convergence  problems  while  using 
the  integration  subroutine  DIVPAG.  This  experience  alerted  us  to  the  possibility  that  the 
venerable  IMSL  subroutine  might  be  the  source  of  the  numerical  noise  we  were  observing. 
Due  to  the  gracious  assistance  of  Prof.  John  Lopez  of  the  Mathematics  Department, 
Pennsylvania  State  University,  we  found  that  DIVPAG  was  not  designed  to  handle  a  NDS 
of  the  type  used  by  our  model;  nor  was  it  reliable  in  working  with  stiff  systems  of 
equations,  which  are  systems  “where  there  are  two  or  more  very  different  scales  of  the 
independent  variable  on  which  the  dependent  variables  are  changing”  (Press  et  al.  1992). 
We  therefore  implemented  the  Numerical  Recipes  semi-implicit  integration  subroutine 
STIFBS  (Press  et  al.  1992),  and  we  were  pleased  that  model  integrations  became 
smoother  and  more  consistent,  not  suffering  eigenvalue  convergence  woes. 

To  attain  the  superior  quality  solutions  afforded  by  STIFBS,  however,  there  was  a 
significant  tradeoff  in  time.  The  new  integration  subroutine  is  considerably  slower  than 
DIVPAG  (nearly  five  times  slower  for  some  cases)  and  this  sluggishness  hindered  research 
progress.  To  alleviate  this,  we  implemented  a  scheme  to  rescale  the  potential  temperature 
perturbation  terms  so  they  more  closely  matched  the  magnitude  of  the  vorticity  terms. 

This  improved  the  balance  of  our  NDS  solutions  by  reducing  the  number  of  iterations 
required  at  each  integration  time  step,  thereby  increasing  the  overall  speed  of  model  runs 
without  affecting  the  linear  stability  results.  To  incorporate  this  time-saving  routine,  we 
redefined  a  new  potential  temperature  perturbation:  ^ 
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^(new)  -  >  (3  -7) 

where  5/ is  the  scale  factor  coefficient.  Instead  of  tediously  multiplying  every  individual 
matrix  in  the  thermodynamic  equation  by  s/,  we  simply  multiplied  the  first  twenty  elements 
of  the  vector  (the  potential  temperature  elements)  by  s/.  To  then  balance  the  equations, 
we  only  had  to  multiply  the  buoyancy  matrix  D  in  the  thermodynamic  equation  by  5/ and 
the  buoyancy  matrix  G  in  the  vorticity  equation  by  (5/)  V  After  integrating  with  these 
alterations,  we  rescale  all  terms  to  their  original  size  to  construct  the  model  solution.  To 
our  pleasure,  the  quality  integrations  given  by  STIFBS  can  now  be  achieved  in  about  half 
the  time. 

We  are  now  ready  to  use  the  model  to  study  very  large  eddy  circulations  of  the 
MABL.  In  Chapter  4,  we  resume  work  on  the  Hi-Res  2  case  study. 
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Chapter  4 
CASE  STUDY 

> 


Since  several  significant  alterations  have  been  made  to  the  model  since  its  last  major 
evaluation  (Lambert  1995),  we  now  ascertain  whether  the  model  solutions  adequately 
represent  observed  atmospheric  phenomena.  To  do  so,  we  compare  and  contrast  our 
numerical  results  with  observations  and  synthetic  aperture  radar  (SAR)  imagery  taken 
during  the  Hi-Res  2  field  study,  which  is  described  in  great  detail  by  Sikora  et  al  (1995) 
and  Lambert  (1995). 

We  perform  the  sequence  of  steps  outlined  in  Figure  2.2  and  Section  2.4  to  arrive  at  a 
numerical  solution  for  the  case  under  investigation.  When  trying  to  capture  via  a  linear 
analysis  the  observed  aspect  ratios  for  an  atmospheric  (as  opposed  to  idealized)  case,  we 
may  use  the  “NO”  loop  in  Figure  2.2  many  times  to  re-select  fi-ee  \rind  parameters  and/or 
the  inversion  strength.  We  iterate  through  this  loop  until  both  of  the  preferred  aspect 
ratios  and  %  match  the  observed  ones  within  +  0. 1  (which  amounts  roughly  to  a  10 
percent  error),  the  specified  error  of  the  SAR  image  spectral  analysis  for  Hi-Res  2  (Shirer 
et  al  1997).  After  a  successful  linear  analysis,  we  determine  the  parameter  values  of  the 
model  using  Hi-Res  2  surface  data,  then  perform  and  analyze  nonlinear  integrations  using 
several  different  scale  factor  values  s/in  (3.7).  This  has  the  effect  of  changing  the  impact 
of  the  nonlinear  dissipation  during  the  integration  without  altering  the  linear  results; 
however,  the  effect  on  the  model  solution  is  substantial,  as  we  demonstrate  in  Section  4.2. 
We  begin  in  Section  4. 1  by  briefly  describing  Hi-Res  2  observations. 
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4. 1 .  Overview  of  Hi-Res  2  Field  Experiment  and  Data 

The  observational  data  used  for  this  case  study  are  from  the  Hi-Res  2  field  experiment 
conducted  in  June  1993  off  the  coast  ofNorth  Carolina,  and  they  are  reported  by  Sikora  et 
al  (1995)  and  Lambert  (1995).  From  16  - 18  June  1993,  two  ships,  the  RA^s  Columbus 
Iselin  and  Bartlett,  took  observations  at  the  sea  surface  and  at  z  =  10  m  within  the  region 
35.2®  -  36.6®  N  latitude  and  73.8®  -  75.0®  W  longitude.  On  17  June  1993,  the  Hi-Res  2 
MABL  measurement  region  coincided  with  that  imaged  by  the  SAR  on  the  European 
Remote  Sensing  Satellite  (ERS-1).  For  this  case  study,  we  choose  a  composite  image  that 
was  taken  at  1538Z  and  was  published  by  Sikora  et  al.  (Figure  1,  1995). 

According  to  Sikora  et  al  (1995),  there  are  two  distinct  patterns  present  on  the  17 
June  1993  SAR  image:  a  marbled  pattern  in  the  northwest  portion  and  a  mottled  pattern 
throughout  the  southeast  portion.  The  transition  between  these  two  regions  occurs  at  the 
Gulf  Stream  North  Wall  (GSNW).  The  marbled  region  north  of  the  GSNW  is  the  SAR 
signature  of  a  stable  MABL,  while  the  mottled  region  south  of  the  GSNW  is  the  signature 
of  a  thermally  unstable  MABL.  Because  we  wish  to  study  secondary  circulations  of  the 
MABL  induced  by  thermal  and/or  inflection  point  instabilities,  we  focus  our  study  on  the 
mottled  region  of  the  SAR  image,  specifically  in  the  vicinity  of  35.8®  N,  74.2®  W,  where 
there  are  discernible  kilometer-scale  alternating  dark/light  patterns.  About  the  time  the 
SAR  image  was  taken,  the  R/V  Columbus  Iselin  was  in  the  area  and  recorded  surface 
observations  that  are  summarized  in  Table  3.1  of  Lambert  (1995).  As  depicted  by  Table 

4.1,  we  regroup  the  17  June  1993  data  into  two  sets.  Sonde  1  and  Sonde  2,  based  upon 
the  times  in  which  the  ship  was  near  the  location  of  each  launch.  Sonde  1,  which  was 
launched  at  1517Z,  represents  the  time-averaged  data  from  1437Z  -  1602Z,  while  Sonde 
2,  which  was  launched  at  1955Z,  represents  the  time  average  from  1941Z  -  2241Z. 
Unfortmately,  boundary  layer  sounding  winds  were  unavailable  for  either  data  set; 
however  R/V  Columbus  Iselin  researchers  reported  that  winds  were  relatively  light 


Table  4.1.  Hi-Res  2  observational  data.  Sonde  1  denotes  data  averaged  between  1437Z  and  1602Z,  and  Sonde  2  denotes  data  averaged 
between  1941Z  and  224 IZ;  both  are  from  17  June  1993.  For  our  study,  we  use  the  Sonde  1  case.  (IW  Columbus  Iselin  data 
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throughout  the  boundary  layer.  Here,  we  use  the  Sonde  1  case  for  our  study  because  it 
was  closer  to  the  time  of  the  SAR  image,  and  we  propose  that  future  work  include  an 
analysis  of  the  Sonde  2  case. 

The  surface  analysis  shown  in  Figure  2  of  Sikora  et  al.  (1995)  depicts  the  synoptic 
conditions  at  1500Z  on  17  June  1993.  Aside  from  a  mesofront  near  the  GSNW  (noted  by 
the  crew  of  the  RA^  Columbus  Iseliri),  the  meteorological  conditions  in  the  region  of 
interest  were  relatively  quiescent:  only  light  winds  and  boundary  layer  cumulus  were 
reported.  This  lack  of  appreciable  synoptic-scale  features  suggests  that  boundary  layer 
convective  eddies  (BLCE)  were  responsible  for  the  stress  patterns  observed  on  the  1538Z 
SAR  image  in  the  region  of  interest  for  this  study  (Sikora  et  al.  1995). 

The  current-relative  winds  at  10  m  for  the  Sonde  1  case  were  0.93  m/s  from  3 19**  and 
the  air-sea  temperature  difference  was  -2.9  **€  (i.e.,  sea  warmer  than  air).  When  plotted 
on  a  Woodcock  (1975)  diagram  (Point  A,  Figure  4.1),  the  data  fall  into  the  three- 
dimensional  regime.  We  therefore  expect  our  model  solution  for  this  particular  case  to  be 
cellular. 

4.2.  Case  Study  Analysis 

We  begin  our  investigation  by  determining  the  values  of  the  lower  boundary  Schramm 
forcing  parameters,  Sj  and  (Figure  2.2,  Step  B)  for  the  Sonde  1  data.  Using  the 
subroutine  described  in  Appendix  A  of  Zuccarello  (1994),  we  input  the  physical  parameter 
values  from  Table  4.1  to  obtain  st  =  1.95  and  s„  =  0.67. 

We  then  conduct  an  eigensystem  linear  analysis  to  determine  the  preferred  values  of 
the  aspect  ratios,  which  yield  the  shape  of  the  horizontal  domain  by  (2.26)  and  (2.27).  For 
the  given  input  parameter  values,  we  iterate  through  a  range  of  relevant  a,  and  ay  values. 


Woodcock  (1975) 


Water  temperature  minus  air  temperature  (C) 


Figure  4. 1.  The  Woodcock  (1975)  diagram  predicts  secondary  circulation  regimes  from  the  current- 

relative  wind  speed  at  10  m  (z*  =  0)  and  the  air/sea  temperature  difference.  Point  A  represents 
the  input  observations  from  the  Hi-Res  2  case  study  (Sonde  1). 
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then  contour  the  real  parts  of  the  50th  eigenvalue  (because  the  eigenvalues  are  sorted  in 
ascending  order  by  the  value  of  their  real  parts).  The  aspect  ratio  region  in  which  the 
largest  positive  real  parts  (or  growth  rates)  are  observed  determines  the  preferred  aspect 
ratios  (Laufersweiler  1987).  While  invoking  the  “NO”  loop  in  Figure  2.2,  we  adjust  the 
free  wind  parameter  values,  inversion  strength  and/or  surface  layer  depth,  and  re-iterate 
through  the  appropriate  Ox  and  ay  values  until  we  approach  the  values  of  the  observed 
aspect  ratios.  We  found  during  our  study  that  this  may  not  be  possible  for  all  cases  in  a 
reasonable  amount  of  time.  Fortunately  for  this  present  study,  after  increasing  from 
0.4  ®C  (Table  4. 1)  to  1.8  ”C,  we  were  able  to  obtain  acceptable  aspect  ratio  values;  the 
linear  stability  analysis  plot  in  Figure  4.2  yields  preferred  aspect  ratios,  ax  and  ay,  that  are 
equal  to  0.7  and  0.8,  respectively.  Shirer  et  a/.(1997)  conducted  a  spectral  analysis  of  the 
1538Z  SAR  image  and  found  ax  =  0.7  +  0.1  and  =  0.9  +  0.1,  as  reported  in  Table  4.1. 
Because  our  findings  are  within  the  error  tolerances  of  the  spectral  analysis,  we  may 
conclude  that  our  NDS  is  capable  of  capturing  the  preferred  horizontal  wavelengths  of  the 
convective  structures  that  are  depicted  on  the  1538Z  SAR  image. 

Table  4.2  summarizes  all  input  parameter  values  used  for  this  Hi-Res  2  case  study. 
Note  that  Re  is  calculated  from  (2.28)  instead  of  being  chosen  directly,  and  that  the 
Coriolis  parameter  f  equals  zero  because  we  expect  the  BLCE  solutions  to  have  time 
scales  on  the  order  of  an  hour  or  less  (Lambert  1995).  Additionally,  we  specify  only  one 
wavenumber  in  each  of  the  x-  and^'-directions  in  the  model,  enabling  us  to  set  /  =  1  and 
/n=l  in  (2.43) -(2.46). 

After  initializing  the  model,  we  integrate  until  we  attain  a  state  of  quasi-equilibration. 
We  observe  this  readily  by  plotting  the  dimensionless  energy  (y^  Ay)  of  the  time- 
dependent  amplitude  coefficients,  which  follows  from  (2.55). 


DETERMINATION  OF  THE 
PREFERRED  ASPECT  RATIO 


Figure  4.2.  Contour  plot  of  the  real  part  of  the  50th  eigenvalue  from  the  linear  stability  analysis  from 
Sonde  1  data.  The  mayimnm  dimensionless  growth  rate  indicates  the  preferred  aspect  ratios 
(0.7, 0.8).  The  associated  model  parameter  values  are  given  in  Table  4.2. 


Table  4.2.  Parameter  values  used  in  case  study  at  point  A  on  Figure  4.1.  36 


PARAMETER  VALUE 


HEIGHT  OF  TOP  OF  DOMAIN  (zr) _ 609 


HEIGHT  OF  BOTTOM  OF  DOMAIN  (/?z^) _ 1^ 


HEIGHT  OF  INVERSION  (zQ _ 409  m 


HEIGHT  OF  SURFACE  LAYER  {zsd _ 2^ 


ASPECT  RATIO  IN  THE  X-DIRECTION  (aj  0.7 


ASPECT  RATIO  IN  THE  Y-DIRECTION  (dr^)  0  g 


REYNOLDS  NUMBER  (Re) _ 359.4 


MEAN  WIND  SPEED  /  DIRECTION  ATz  =  Zr _ 6.0  ms^  /  300° 


MEAN  WIND  SPEED  /  DIRECTION  AT  z  = _ 0.93 /W5~V  319° 


U  MEAN  WIND  PROFILE  PARAMETER  (du) _ ^ 


V  MEAN  WIND  PROFILE  PARAMETER  (dy) _ -3.0 


ACCELERATION  OF  GRAVITY  (g) _ 9.8  ms'^ 


CORIOLIS  PARAMETER  (/) _ 0^ 


EDDY  VISCOSITY  (v)  1 


EDDY  CONDUCTIVITY  (x:)  10 


PRANDTL  NUMBER  (F)  a  1 


SEA  SURFACE  TEMPERATURE  (Too) _ 27.8  °C 


TEMPERATURE  AT  z^Hlb  (Tip) _ _ 24.9 


MIXING  RATIO  AT  SEA  SURFACE 


MIXING  RATIO  ATz  =  /iib 


RAYLEIGH  NUMBER  (Fa) _ _ ^ 


SCHRAMM  VELOCITY  CONSTANT  (s„) _ 0.67 


SCHRAMM  THERMAL  CONSTANT  (sr) 


1.95 
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As  noted  in  Section  3. 1,  we  empirically  chose  the  value  of  the  numerical  nonlinear 
dissipation  coefficient  to  be  0.1;  it  remains  fixed-as  do  the  parameters  in  Table  4.2-for 
the  two  Sonde  1  cases  presented  here.  What  we  do  alter,  however,  is  the  value  of  the 
scale  factor  j/from  (3.7),  which  was  initially  implemented  to  expedite  the  integration 
process  in  a  way  that  has  absolutely  no  impact  on  the  linear  analysis,  nor  in  principle,  on 
the  nonlinear  advection  terms.  However,  changing  the  value  of  j/(and  nothing  else)  does 
alter  the  weight  of  the  nonlinear  dissipation  in  the  integration  to  yield  rather  interesting 
results.  Table  4.3  shows  the  two  Sonde  1  cases  that  we  integrated  with  the  scale  factor 
and  run  time  shown. 


Table  4.3.  Two  Sonde  1  cases  determined  by  scale  factor  value. 


Case 

d 

Sf  _ 

Length  of  Integration 

1 

0.1 

10’^ 

175,000  s 

2 

0.1 

10'^ 

150,000  s 

We  arrived  at  these  values  after  experimenting  with  d  values  in  the  range  of  0  to  10, 
and  .^values  in  the  range  of  10“*  to  1 .  For  the  case  with  no  dissipation  (d=0,Sf=  1),  the 
dimensionless  energy  began  to  blow  up  at  175,000  s;  in  contrast,  cases  with  strong 
numerical  dissipation  displayed  extremely  long,  periodic  oscillations  with  no  sign  of  quasi¬ 
equilibration  within  reasonable  periods  of  time.  For  each  order  of  magnitude  decrease  in  Sf 
(e.g.,  1  to  10’*),  we  realized  an  approximate  10  percent  decrease  in  integration  time.  The 
selections  of  d  and  5/ values  shown  in  Table  4.3  represent  what  we  found  to  be  the  best 
combinations  for  speed  and  run  quality. 

The  dimensionless  energy  is  plotted  for  Case  1  and  Case  2  in  Figures  4.3  and  4.4, 
respectively.  Case  2  reaches  a  state  of  quasi-equilibration  around  the  dimensionless 
energy  value  7.4  x  10^**,  and  so  we  stopped  the  run  at  150,000  s  (approximately  42  hours). 


DIMENSIONLESS  ENERGY 


2E+11 


1.6E+11 


1.2E+11 


8E+10 


4E+10 


DIMENSIONLESS  ENERGY 
(CASE  1) 


8000  12000 

TIME  (s)x10 


Figure  4.3.  Case  1  time  series  plot  of  the  dimensionless  energy  for  the  50  time-dependent  amplitude 
coefBcients.  The  system  is  quasi-equilibrated. 
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In  contrast,  Case  1  appears  to  be  settling  near  a  value  of  9.6  x  10^°  by  150,000  s. 

Although  the  oscillation  amplitude  increases  in  magnitude  between  150,000  s  and  175,000 
s,  the  mean  energy  is  still  less  than  1.0  x  lO”.  We  integrate  Case  1  further  to  200,000  s  to 
confirm  that  the  oscillations  do  not  continue  to  intensify  without  bound,  and  as  Figure  4.3 
shows,  they  apparently  do  not.  They  initially  increase  in  magnitude,  then  decrease,  with 
oscillations  centered  about  1.1  x  10“.  We  choose  175,000  s  as  the  case  study  termination 
point,  however,  because  the  model  solution  appears  to  stabilize  by  then,  with  minimal 
transient  activity.  Time  series  plots  of  the  potential  temperature  amplitude  coefficients  9^^ 
for  each  case  (Figure  4.5  for  Case  1;  Figure  4.6  for  Case  2)  show  regular  periodic 
oscillations,  indicating  that  our  choices  for  run  termination  times  are  appropriate.  Note 
that  Case  1  and  Case  2  have  dominant  solution  periods  of  15  minutes  and  67  minutes, 
respectively,  which  are  consistent  with  the  findings  of  Shirer  (1986),  Pyle  (1987),  and 
Haack  and  Shirer  (1992)  for  BLCE. 

Because  the  ultimate  goal  of  this  project  is  to  study  the  impact  that  MABL  secondary 
circulations  have  on  the  sea  surface  stress  patterns  as  observed  by  SAR,  we  calculate  the 
stress  at  the  lower  boundary  of  the  model  domain  using  the  following  equation: 

r,  =  Qft[(C/.„+»,)'+(F.„+v,^)].  (4.1) 

Here,  and  are  the  horizontal  components  of  the  mean  background  wind  at  the 
lower  boundary;  and  v,.  are  the  horizontal  components  of  the  wind  due  to  the  roll  or 
cell  at  the  lower  boundary  /jlb;  is  the  density  from  (2.2)  with  a  typical  value  of  1  kg/m^-y 
and  Q  is  the  coefficient  of  drag,  with  a  typical  value  of  1/900  (Stull  1988).  The 
computational  output  (CO)  subroutine  of  the  model  also  calculates  the  stress  induced 
strictly  by  and  v,.  to  display  the  contribution  from  the  rolls  or  cells  independent  of  the 
mean  wind  contribution.  For  the  stress  plots  presented  in  this  study,  we  chose  to  include 
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POTENTIAL  TEMPERATURE  AMPLITUDE 
COEFFICIENT  (CASE1) 


Figure  4.5.  Time  series  (last  15,000  s)  of  one  of  the  time-dependent  potential  temperature  amplitude 
coefficients  ( ^, , )  for  Case  1,  demonstrating  that  the  solution  has  a  dominant  period  of 
approximately  15  minutes. 


TEMP  COEFFICIENT  VALUE 
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POTENTIAL  TEMPERATURE  AMPLITUDE 
COEFFICIENT  (CASE  2) 


Figure  4.6.  Time  series  (last  20,000  s)  of  one  of  the  time-dependent  potential  temperature  amplitude 
coefficients  ( )  for  Case  2,  demonstrating  that  the  solution  has  a  dominant  period  of 
approximately  67  minutes. 
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the  contribution  of  the  background  wind,  as  that  would  apply  better  to  the  atmospheric 
case. 

A  salient  application  of  stress  plan\dews  is  the  determination  of  the  structural  type  of 
BLCE  regime  under  investigation,  whether  quasi-two-dimensional  with  rolls,  fully  three- 
dimensional  with  cells,  or  somewhere  in  between.  While  we  are  primarily  concerned  with 
the  regime  of  the  final  stabilized  solution  for  each  case,  we  are  also  interested  in  observing 
the  pattern  evolution  between  regimes.  To  do  so,  we  occasionally  extract  data  during  the 
course  of  the  integrations  and  pass  it  through  the  CO  subroutine.  On  the  dimensionless 
energy  time  series  plots  (Figures  4.3  and  4.4),  we  have  identified  the  points  at  which  data 
were  extracted~A  through  E  for  Case  1,  and  W  through  Z  for  Case  2-and  we  have 
summarized  the  results  in  Table  4.4.  As  we  illustrate  below,  note  that  in  Case  1  the  model 


Table  4.4.  Summary  of  structure  types  at  the  points  depicted  on  the  energy  plots  (Figures  4.3  and  4.4) 
for  Case  1  and  Case  2.  The  type  is  determined  from  the  surface  stress  patterns. 


CASE 

POINT 

STRUCTURE 

ENERGY  LEVEL 

1 

25,000 

A 

Rolls 

1.41x10“ 

1 

75,000 

B 

Transition 

1.59x10“ 

1 

100,000 

C 

CeUs 

1.48x10“ 

1 

150,000 

D 

Cells 

9.58  X  10‘“ 

1 

175,000 

E 

CeUs 

1.11x10“ 

2 

75,000 

W 

RoUs 

7.65  X  10‘® 

2 

100,000 

X 

RoUs 

6.75  X  10‘® 

2 

125,000 

Y 

Rolls 

6.86x10^" 

2 

150,000 

Z 

RoUs 

7.44x10*” 

solution  starts  as  quasi-two-dimensional,  passes  through  a  transitional  phase  and  then 
becomes  fiilly  three-dimensional.  Case  2,  in  contrast,  is  strictly  quasi-two-dimensional 
throughout  the  length  of  the  integration,  with  no  indications  of  changing.  Remarkably, 
both  cases  were  run  fi’om  the  same  initial  conditions.  Because  the  magnitude  of  the 
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dimensionless  energy  is  much  greater  overall  in  Case  1,  we  speculate  that  the  transition 
from  quasi-two-dimensional  to  three-dimensional  patterns  requires  higher  levels  of 
energy — energy  that  is  reduced  in  Case  2  by  using  a  larger  value  of  Sf,  which  serves  to 
strengthen  the  nonlinear  dissipation.  A  majority  of  the  preliminary  runs  we  performed 
produced  energy  signatures  like  that  depicted  in  Case  1  and  yielded  three-dimensional 
solutions  after  quasi-equilibration  had  occurred.  We  therefore  postulate  that  higher 
energy  levels  are  a  necessary  (but  not  sufficient)  condition  for  a  regime  to  transition  from 
a  quasi-two-  to  three-dimensional  pattern.  Owing  to  time  limitations,  we  were  unable  to 
explore  this  hypothesis  fiirther  but  suggest  that  future  studies  address  it. 

The  plot  of  stress  magnitudes  at  the  lower  boundary  for  Case  1  (at  Point  E,  Figure 
4.3)  is  shown  in  Figure  4.7.  Higher  stress  regions  are  depicted  by  lighter  shading  to  match 
the  brightness  of  the  high-stress  regions  on  a  SAR  image,  and  we  are  gratified  that  the 
magnitudes  of  the  dimensional  stress,  with  orders  of  magnitude  of  10'^  Nm'^,  are 
consistent  with  those  found  by  Alpers  and  Brummer  (1994)  (although  their  10  m  winds 
were  more  intense).  As  previously  noted,  the  Case  1  stress  pattern  clearly  exhibits  a  fully 
three-dimensional  flow  field  whose  structure  is  oriented  nearly  perpendicular  to  the 
direction  of  the  mean  background  wind  at  the  lower  boundary.  We  contrast  this  result 
with  that  in  Figure  4.8,  which  shows  the  surface  stress  for  Case  1  at  100,000  s  (Point  C, 
Figure  4.3).  Here  the  orientation  is  aligned  14®  to  the  right  of  the  background  wind, 
which  is  well  within  the  range  of  typically  observed  values  of  -20®  to  +30®  (Etling  and 
BroAvn  1993).  This  orientation  also  more  closely  matches  the  pattern  depicted  by  the 
1538Z  SAR  image  shown  in  Sikora  etal  (1995)  as  well  as  Lambert’s  (1995)  findings. 
From  a  modeling  standpoint,  discovering  that  cells  can  become  oriented  normal  to  the 
background  wind  is  a  success;  this  is  the  first  case  in  which  we  observed  such  behavior, 
assuring  us  that  the  orientation  is  not  locked  to  the  direction  of  the  mean  wind  as  a  result 
of  a  modeling  assumption.  From  a  meteorological  perspective,  however,  we  felt  initial 
disappointment  because  we  had  not  captured  the  same  cell  orientation  as  the  SAR  image 
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STRESS  AT  LOWER  BOUNDARY 
(CASE  1 ,  POINT  E) 


Figure  4.7.  Planview  of  stress  magnitudes  at  z*  =  0  imposed  by  the  background  wind  and  secondary 

circulations  (cells),  as  predicted  in  Figure  4.1.  Consistent  with  SAR  imagery,  regions  of  higher 
stress  are  shaded  lightly  and  regions  of  lower  stress  are  shaded  heavily.  The  domain  shape 
(x*/y*)  is  chosen  to  agree  with  the  ratio  ajox  =  1.15,  and  Xr*  and  Yr*  represent  the  rotated  cell 
coordinate  system.  The  along-mean  wind  axis  Xr*  is  per^ndicular  to  the  cross-mean  axis  Yr* 
but  does  not  appear  as  such  due  to  the  domain  shape  (x  /y  )  scaling.  The  contour  interval  is 
5  X  10"^  and  the  +  sign  denotes  the  point  at  which  profiles  are  calculated. 
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STRESS  AT  LOWER  BOUNDARY 
(CASE  1 ,  POINT  C) 


Figure  4.8.  Planview  of  stress  magnitudes  at  z*  =  0  imposed  by  the  background  wind  and  secondary 
circulations  (cells),  as  predicted  by  Figure  4.1.  Consistent  with  SAR  imagery,  regions  of  high 
stress  are  shaded  lightly  and  regions  of  lower  stress  are  shaded  heavily.  The  domain  shape 
(xVy*)  is  chosen  to  agree  with  the  ratio  a/a,  =  1.15,  and  the  contour  interval  is  the  same  as  that 
used  in  Figure  4.7. 
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at  the  same  time  that  model  transients  had  settled  down.  However,  a  SAR  image  records 
an  instant  in  time,  not  an  evolution,  and  it  is  quite  possible  that  a  change  occurred  in  the 
orientation  of  the  sea  surface  stress  pattern  after  the  1538Z  SAR  image  was  taken.  Of 
course  we  cannot  prove  such  an  assertion  here,  but  we  note  that  having  the  secondary 
circulation  aligned  normal  to  the  direction  of  the  mean  wind  is  not  without  precedent  for 
an  inflection  point  mode  BLCE  (e.g..  Brown  1980).  Later,  we  will  show  that  the  Case  2 
modified  (post-cell)  wind  speed  profile  contains  two  points  of  inflection,  consistent  with 
our  findings  here. 

In  this  study,  we  do  not  focus  on  the  transition  from  quasi-two-  to  three-dimensional 
regimes;  however,  we  include  in  Figure  4.9  the  stress  pattern  plot  for  Point  B  of  Figure 
4.3,  which  illustrates  the  transitional  phase  stress  pattern  for  Case  1.  Note  that  while  the 
quasi-two-dimensional  roll  structure  is  still  quite  prominent,  there  is  evidence  of  cellular 
activity. 

The  stress  pattern  for  Case  2  (Point  Z,  Figure  4.4)  is  depicted  in  Figure  4. 10  and  is 
quasi-two-dimensional,  indicating  that  rolls  are  dominant  in  the  model  solution.  There  is  a 
hint  of  three-dimensionality  in  the  stress  pattern  (small  cellular  structures  along  the 
diagonal),  but  it  is  negligible;  the  phenomenon  may  be  due  to  limitations  of  the  contouring 
software  in  handling  the  very  small  variations  along  the  SE-NW  diagonal,  or  simply  that 
the  gridding  used  was  too  coarse.  The  rolls  for  this  case  are  closely  aligned  with  the 
direction  of  the  mean  background  wind  and  remain  fixed  in  orientation  throughout  the 
length  of  the  run,  indicating  a  more  robust  model  solution.  Based  upon  the  extensive 
work  done  by  Woodcock  (1940,  1975)  and  DeardorfiF  (1976),  we  recall  fi’om  Figure  4. 1 
that  the  atmospheric  parameter  values  used  (Table  4. 1)  should  yield  a  three-dimensional 
cellular  solution.  Accordingly,  we  attribute  the  quasi-two-dimensionality  observed  in  Case 
2  strictly  to  numerics—the  larger  j/ value  simply  altered  dissipation  sufficiently  to  cause  the 
model  solution  to  seek  a  different,  lower  energy  attractor,  a  quasi-two-dimensional  one. 
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STRESS  AT  LOWER  BOUNDARY 
(CASE  1,  POINT  B) 


0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0 


X* 


Figure  4.9.  Planview  of  stress  magnitudes  at  z*  =  0  imposed  by  the  background  wind  and  secondary 

circulations  in  transition  between  rolls  and  cells.  Consistent  with  SAR  imagery,  regions  of  high 
stress  are  shaded  lightly  and  regions  of  lower  stress  are  shaded  heavily.  The  domain  shape 
(x7y*)  is  chosen  to  agree  with  the  ratio  Oy/Ox  =  1. 15,  and  the  contour  interval  is  the  same  as  that 
used  in  Figure  4.7. 


STRESS  AT  LOWER  BOUNDARY 
(CASE  2,  POINT  Z) 
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Figure  4. 10.  Planview  of  stress  magnitudes  at  z*  =  0  imposed  by  the  background  wind  and  secondary 

circulations  (rolls).  Consistent  with  SAR  imagery,  regions  of  hi^  stress  are  shaded  lightly  and 
regions  of  lower  stress  are  shaded  heavily.  The  domain  shape  (x/y')  is  chosen  to  agree  with  the 
ratio  Oj/Ox  =  1.15,  and  Xr*  and  Yr*  represent  the  rotated  roll  coordinate  system  as  in  Figure  4.7. 
The  contour  interval  is  also  the  same  as  that  used  in  Figure  4.7,  and  the  +  sign  denotes  the  point 
at  which  profiles  are  calculated. 
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We  therefore  postulate  that  Case  1  is  more  physically  plausible  than  Case  2  with  the 
atmospheric  input  variables  used.  Nevertheless,  we  will  compare  and  contrast  both  cases 
to  further  evaluate  model  capabilities  and  limitations. 

Because  the  cellular  structures  in  Case  1  are  aligned  nearly  normal  to  the  direction  of 
the  background  mean  wind  (Figure  4.7)  and  because  the  rolls  in  Case  2  are  nearly  parallel 
to  the  mean  wind,  we  expect  the  dominant  solution  period  to  be  smaller  for  Case  1. 
Agreeably,  it  is,  as  we  determined  from  the  time  series  plots  for  (Figures  4.5  and  4.6); 
the  Case  1  period  is  15  minutes  while  the  Case  2  period  is  67  minutes. 

On  Figures  4.7  and  4.9  we  have  drawn  new,  rotated  axes,  and  Yr*,  along  which 
we  take  cross  sections  of  the  domain  to  contour  values  of  the  streamfunctions,  vertical 
velocity  and  potential  temperature.  These  rotated  axes  are  determined  differently  from  the 
X*  and  3?* axes  used  in  Lambert  (1995).  Here,  Xr*  and  Yr*  represent  the  cross-  and 
along-mean  wind  axes,  and  they  are  perpendicular  to  one  another.  On  the  stress 
planviews,  however,  they  do  not  appear  as  such  due  to  domain  shape  scaling.  The  rotated 
axes  are  given  by 


(4.2) 

(4.3) 

and  the  rotated  streamfunctions  are  specified  by 


Xr*=^x*  -y* 
Yr*  =  x*  , 


«  ♦ 
=  v  -1  ■ 


(4.4) 

(4.5) 
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Note  that  we  only  illustrate  physically  relevant  streamfunction  components  in  our  cross 
sections;  that  is,  7*  is  contoured  only  along  the  Xr*  axis  as  7* ,  and  y/  is  contoured  only 
along  the  Yr*  axis  as  y/* .  For  Case  1,  Xr*  runs  along  the  orientation  of  the  cells  and  Yr* 
runs  across  the  cells.  For  Case  2,  Xr*  and  Yr*  are  the  cross-  and  along-roll  axes, 
respectively. 

We  begin  by  contouring  values  of  the  streamfunctions  for  both  cases  (Figures  4.11  and 
4.12)  to  ensure  dimensional  consistency  with  the  stress  plots  and  to  determine  the  vertical 
extent  of  the  secondary  circulations.  A  cursory  glance  at  the  flows  shown  for  Case  1 
(Figure  4. 1 1)  clearly  suggests  a  fully  three-dimensional  structure  because  flows  are  evident 
in  both  cross  sections— this  is  what  we  expected  firom  Figure  4.7.  Similarly,  the  quasi-two- 
dimensional  structure  displayed  by  the  stress  pattern  in  Figure  4.10  for  Case  2  clearly 
coincides  with  the  flows  visualized  in  Figure  4. 12,  for  there  are  only  circulations  present 
along  the  Ar*  axis,  the  cross-roll  axis  in  this  case.  Overall,  the  circulations  are  stronger  in 
Case  1  than  in  Case  2. 

We  had  hoped  that  the  addition  of  an  inversion  (at  z*  =  0.7)  would  cap  the  circulations 
as  it  did  for  Laufersweiler  and  Shirer  (1995).  From  Figures  4. 1 1  and  4. 12,  we  can  see  that 
it  appears  to  be  doing  so  for  both  cases,  with  only  slight  anomalies.  For  Case  1  (Figure 
4. 1 1),  the  along-Ar*  axis  plot  shows  circulations  capped  almost  perfectly  by  the  inversion 
with  very  weak,  and  certainly  acceptable,  perturbations  above  the  inversion  base  at  z  = 
0.7.  Along  the  Yr*  section,  the  circulations  penetrate  into  the  inversion,  though  the 
relative  weakness  of  y/* above z*  is  comforting.  Given  the  relatively  small  inversion 
strength  (1.8  ®C/200m),  this  too  seems  physically  acceptable.  Similarly,  the  Case  2  flow 
shown  along  the  Xr*  axis  in  Figure  4. 12  is  quite  reasonable  with  a  well-defined  circulation 
below  z*  and  only  negligible  perturbations  into  the  inversion.  Overall,  these  cross  sections 
are  encouraging  and  illustrate  that  the  model  is  capable  of  adequately  capping  weak  to 
moderate  flows  with  an  inversion. 


STREAMFUNCTION 
CROSS  SECTIONS 
(CASE  1) 


Figure  4. 11.  Cross  sections  of  contoured  streamfimctions  using  rotated  horizontal  axes  from  Figure  4.7; 

*  * 

(Top)  T}  is  contoured  along  the  axis,  and  (Bottom)  is  contoured  along  the  Yr*  axis. 

Shaded  regions  represent  negative  values  while  lighter  regions  represent  positive  values.  The 
contour  interval  is  consistent  in  the  top  and  bottom  cross  sections,  as  well  as  with  Figure  4. 12. 


STREAMFUNCTION 
CROSS  SECTIONS 
(CASE  2) 
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Figure  4. 12.  Cross  sections  of  contoured  streamfunctions  using  rotated  horizontal  axes  from  Figure  4.9; 

(Top)  is  contoured  along  theXr*  axis,  and  (Bottom)  is  contoured  along  the  Yr*  axis. 

Shaded  regions  represent  negative  values  while  lighter  regions  represent  positive  values.  The 
contour  interval  is  consistent  in  the  top  and  bottom  cross  sections,  as  well  as  with  Figure  4.11. 
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Figures  4. 13  and  4. 14  display  cross  sections  of  vertical  velocities  for  Cases  1  and  2, 
respectively.  One  notable  feature  in  both  cases  is  the  presence  of  vertical  velocity  maxima 
at  approximately  0.3  of  the  depth  of  the  boundary  layer  (z  «  0.2),  which  is  consistent 
with  convective  boundary  layer  observations  (e.g.,  Lenschow  et  al.  1980)  and  other  more 
complex  model  simulations  (e.g.,  large  eddy  simulation  [LES],  Schmidt  and  Schumann 
[1989]).  In  the  Yr*  cross  section  of  Case  1  (Figure  4. 13)  and  both  cuts  from  Case  2 
(Figure  4.14),  the  vertical  velocity  contours  are  contained  below  the  base  of  the  inversion 
at  Zi  -  0.7,  with  the  exception  of  very  minor  perturbations.  The  Xr*  cross  section  of  Case 
1,  however,  appears  as  an  anomaly  with  a  wavenumber-three  signature  and  the  maximum 
updraft  just  below  r,.  We  attribute  this  pattern  to  the  location  of  the  Xr*  axis  in  Figure 
4.7,  which  cuts  through  a  complicated  three-dimensional  pattern  between  the  high-  and 
low-stress  regions,  and  encounters  a  plethora  of  small-scale  contour  variations.  Had  we 
taken  a  different  cross  section,  we  would  expect  cleaner  results.  We  are  quite  pleased 
with  the  Xr*  cross  section  of  Case  2,  for  it  clearly  captures  the  observed  phenomenon  of 
narrow  updrafts  and  broad  downdrafts.  We  are  somewhat  perplexed  by  the  Yr*  section  of 
Case  2,  however,  since  the  existence  of  updrafts  in  this  cross  sections  suggests  a  three- 
dimensional  flow;  yet  all  other  indicators  point  to  a  quasi-two-dimensional  regime.  We 
defer  speculation  until  the  solution  is  visualized  in  three  dimensions. 

We  turn  now  to  vertical  velocity  profiles  for  a  snapshot  look  at  the  variation  of  w  with 
z*  at  a  fixed  horizontal  point  (x*  =  0.8,  y*  =  0.5),  denoted  on  Figures  4.7  and  4. 10.  A  Case 
1  profile  drawn  in  Figure  4.15  depicts  a  maximum  updraft  just  above  z*  =  0.2  (consistent 
mth  the  cross  section  results),  a  downdraft  just  above  z*  =  0.5  and  a  secondary,  weaker 
updraft  in  the  inversion  layer.  We  believe  the  two  former  attributes  accurately  represent 
the  physics  of  the  boundary  layer,  while  the  latter  attribute  shows  evidence  of  a  very  slight 
circulation  within  the  inversion  layer,  consistent  with  the  streamfimction  and  vertical 
velocity  perturbations  noted  in  Figures  4. 1 1  and  4. 13,  respectively.  The  Case  2  vertical 
velocity  profile  is  shown  in  Figure  4. 16  and  represents  a  nearly  ideal  profile,  for  there  is 


VERTICAL  VELOCITY 
CROSS  SECTIONS 
(CASE  1) 
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Figure  4.13.  Cross  sections  of  contoured  vertical  velocity  fields  using  rotated  horizontal  axes  from  Figure 
4.7:  (Top)  w*  is  contoured  along  the  Xr*  axis,  and  (Bottom)  w*  is  contomed  along  the  Yr*  axis. 
Shaded  regions  represent  negative  values  (downdrafts)  while  lighter  regions  represent  positive 
values  (updrafts).  The  contour  interval  is  consistent  in  the  top  and  bottom  cross  sections,  as  well 
as  with  Figure  4. 14. 
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VERTICAL  VELOCITY 
CROSS  SECTIONS 
(CASE  2) 


Figure  4. 14.  Cross  sections  of  contoured  vertical  velocity  fields  using  rotated  horizontal  axes  from  Figure 
4.9:  (Top)  w*  is  contoured  along  the  Yr*  axis,  and  (Bottom)  w*  is  contoured  along  the  Yr*  axis. 
Shaded  regions  represent  negative  values  (downdrafts)  while  lighter  regions  represent  positive 
values  (updrafts).  The  contour  interval  is  consistent  in  the  top  and  bottom  ctoss  sections,  as  well 
as  with  Figure  4.13. 
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only  one  maximum  in  the  vertical  column  at  approximately  0.3  times  the  depth  of  the 
boundary  layer.  Additionally,  w  becomes  negligible  just  below  the  base  of  the  inversion 
and  remains  small  throughout  the  inversion  layer.  We  are  also  encouraged  that  the 
magnitudes  of  w  found  for  both  cases  are  acceptable  for  BLCE  (Etling  and  Brown  1993). 

We  now  examine  the  potential  temperature  perturbation  cross  sections  shown  in 
Figures  4.17  and  4. 18  for  Cases  1  and  2,  respectively.  A  cursory  glance  reveals  that  the 
model  has  difficulty  resolving  temperature  in  these  cross  sectional  views.  However,  there 
are  three  noteworthy  phenomena  captured  by  the  model;  1)  The  region  where  the  greatest 
temperature  perturbations  occur  is  at  the  lower  boundary,  indicating  that  our  boundary 
conditions  have  captured  the  physics  of  a  warm  ocean  below  a  cooler  MABL,  as  Schmidt 
and  Schumann  (1989)  found  with  LES.  2)  There  is  a  clear  demarcation  in  potential 
temperature  between  the  boundary  layer  and  the  inversion  layer,  and  we  take  comfort  that 
the  line  where  del’ jc!:*  =0  is  in  close  proximity  to  z*  =  0.7.  3)  The  dimensional  values  of 
the  potential  temperature  perturbations  range  from  order  0.1  AT  to  1.0  AT;  the  lower  value  is 
in  accordance  with  the  results  of  Haack  and  Shirer  (1992). 

The  wavenumber  3  and  4  patterns  observed  have  only  a  minor  effect  on  the  vertical 
transport  of  thermal  energy  in  the  model  (Figures  4.19  and  4.20).  These  convincing 
dimensionless  heat  flux  profiles  reveal  small  negative  heat  fluxes  in  the  inversion  layer,  as 
is  typically  observed.  Both  figures  also  show  maximum  fluxes  near  the  bottom  of  the 
domain,  which  is  to  be  expected  when  there  is  heating  from  below.  The  magnitudes  of  the 
heat  fluxes  are  of  order  XQ'^Kms'  for  Case  1,  and  order  Kms'  for  Case  2;  the  smaller 
values  obtained  in  Case  2  are  consistent  with  the  findings  of  Laufersweiler  and  Shirer 
(1995),  who  also  studied  a  quasi-two-dimensional  case.  The  Case  2  heat  flux  profile 
(Figure  4.20)  is  especially  encouraging,  for  it  decreases  nearly  linearly  from  a  relatively 
large  positive  value  at  r*  =  0. 15  to  a  relatively  smaU  negative  value  at  z,*;  then  above  z,*,  it 
increases  in  value  towards  zero,  showing  the  entrainment  (between  z  =  0.6  and  z  =  0.7) 


POTENTIALTEMPERATURE 
CROSS  SECTIONS 
(CASE  1) 
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Figure  4. 17.  Cross  sections  of  contoured  potential  temperature  fields  using  rotated  horizontal  axes  from 
Figure  4.7;  (Top)  9*  is  contoured  along  the  A>*  axis,  and  (Bottom)  9*  is  contoured  along  the  Yr* 
axis.  Shaded  regions  represent  negative  values  while  lighter  regions  represent  positive  values. 
The  contour  interval  is  consistent  in  the  top  and  bottom  frames,  as  well  as  with  Figure  4. 18. 
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Figure  4. 18.  Cross  sections  of  contoured  potential  temperature  fields  using  rotated  horizontal  axes  from 
Figure  4.9:  (Top)  9*  is  contoured  along  the  Xr*  axis,  and  (Bottom)  6*  is  contoured  along  the  Yr* 
axis.  Shaded  regions  represent  negative  values  while  lighter  regions  represent  positive  values. 
The  contour  interval  is  consistent  in  the  top  and  bottom  frames,  as  well  as  with  Figure  4.17. 
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of  warm  inversion  layer  air  into  the  boundary  layer.  At  the  lower  boundary,  however,  the 
fluxes  return  to  near  zero  values,  which  is  not  representative  of  a  cool  MABL  over  a 
warm  ocean  when  the  proper  subgrid  components  are  included.  The  heat  flux  reverts  to 
zero  as  a  consequence  of  the  imposed  boundary  condition  (2.36),  but  is  acceptable  for  the 
BLCE  scale  (e.g.,  Briimmer  1985;  Becker  1987;  Chlond  1987). 

The  vertical  transport  of  horizontal  momentum  (i.e.,  the  mean  momentum  flux)  is 
depicted  in  Figure  4.21  for  Case  1  and  Figure  4.22  for  Case  2.  Though  lacking  the  detail 
given  by  calculating  the  momentum  in  the  separate  x*-  and>'*-components,  these  mean 
momentum  profiles,  which  are  calculated  using  standard  flux  equations  (Stull  1988),  are 
good  indicators  of  the  overall  momentum  transport  by  the  model  solutions.  From  Figure 
4,22,  we  observe  that  the  maximum  flux  of  momentum  occurs  near  z  —  0.35  for  Case  2. 
This  is  reassuring,  for  it  is  relatively  consistent  with  our  findings  for  the  maximum  vertical 
velocity,  discussed  earlier;  we  are  also  pleased  that  the  flux  values  in  the  inversion  are 
negligible.  Case  1,  which  is  shown  in  Figure  4.21,  shows  a  local  maximum  at  z  =  0.3  but 
a  much  greater  maximum  in  the  inversion  layer.  While  this  correlates  well  with  the  Case  1 
vertical  velocity  profile  in  Figure  4.15,  as  we  expect,  the  presence  of  a  momentum  flux 
mayimnm  in  the  inversion  is  somewhat  disconcerting  and  worthy  of  further  study. 

The  secondary  circulations  determined  by  the  model  alter  the  background  mean  wind 
profile  through  horizontally-independent  modifications,  which  come  only  fi’om  the 

and  components  (Haack  and  Shirer  1992).  In  Figure  4.23  (Case  1)  and  Figure  4.24 

(Case  2),  we  plot  the  background  wind  speed  |  V(z)|  and  the  modified  plus  background 
wind  speed  lV(z)  +  Vmod(z)|  on  the  same  graph.  This  latter  wind  value  at  {z  =  hi^,  known 
also  as  the  post  roll/cell  wind,  is  what  would  be  diagnosed  fi'om  a  SAR  image.  For  both 
cases,  we  clearly  show  that  the  BLCE  circulations  are  redistributing  the  shear  in  the 
boundary  layer,  as  was  noted  by  Haack  and  Shirer  (1992).  We  also  note  that  the 
magnitude  of  the  change  in  wind  speed  is  in  line  with  the  findings  of  Haack  and  Shirer 
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Figure  4.23.  Profile  of  background  and  modified  (post-cell)  wind  speeds  for  Case  1 
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(1992),  and  is  sufficient  to  alter  the  linear  analysis  significantly.  The  modified  wind  speed 
profiles  are  quite  similar  for  Cases  1  and  2,  although  the  wind  speed  decreases  more 
rapidly  near  the  surface  for  Case  1,  indicating  slightly  greater  shear  in  the  lowest  portions 
of  the  MABL.  According  to  the  Woodcock  (1975)  diagram  in  Figure  4. 1,  we  would 
expect  regimes  with  stronger  winds  at  10  /w  {hu^  to  have  a  higher  likelihood  of  becoming 
quasi-two-dimensional  than  regimes  with  weaker  lOm  winds.  This  is  consistent  with  our 
case  study  result  because  the  modified  wind  speed  at  hie  for  Case  1  is  approximately  0.4 
m/s,  while  the  speed  at  Hlb  for  Case  2  is  1.3  m/s.  Though  not  a  large  absolute  difference, 
the  relative  difference  is  not  insignificant.  We  therefore  postulate  that  Case  1  evolved  into 
a  three-dimensional  solution  consistent  with  the  reduction  in  the  wind  speed  at  the  lower 
boundary,  while  Case  2  remmned  a  quasi-two-dimensional  solution  consistent  with  an 
increase  in  the  wind  speed  at  the  lower  boundary. 

We  further  propose  that  the  Woodcock  (1975)  diagram  be  reformulated  as  a  statistical 
probability  graph  rather  than  an  absolute  predictor  of  regime  t5q)e.  We  could  do  this  by 
adding  contours  of  constant  probability  to  the  chart  to  determine  the  odds  of  obtaining 
rolls  or  cells  based  upon  the  input  parameter  values  (i.e.,  location  on  the  diagram);  such 
contours  could  be  resolved  experimentally  based  upon  the  relative  number  of  runs  that 
give  quasi-two-  or  three-dimensional  solutions.  For  example,  the  demarcation  line 
between  rolls  and  cells  on  Figure  4. 1  would  most  likely  represent  the  50  percent 
probability  line  of  obtaining  one  or  the  other  regime  type.  Our  early  model  runs  also 
support  a  probabilistic  approach  to  interpretation  of  the  diagram  (e.g..  Table  4.4  for  Case 
1).  In  Figures  4.23  and  4.24,  we  observe  that  shear  is  minimal  in  the  inversion  layer, 
consistent  with  observations.  Additionally,  the  modified  wind  speed  profiles  for  each  case 
contain  two  points  of  inflection,  whereas  the  background  wind  speed  profile  contained 
none.  We  stated  previously  that  an  inflection  point  in  the  wind  speed  profile  is  a  necessary 
condition  for  cells  to  align  normal  to  the  direction  of  the  mean  wind  (Brown  1980),  as  we 
found  for  Case  1;  however,  it  is  not  a  sufficient  condition  because  our  Case  2  rolls  aligned 
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nearly  parallel  to  the  direction  of  the  mean  wind  despite  the  presence  of  inflection  points  in 
the  modified  wind  profile. 

The  wind  direction  profiles,  though  not  shown  here,  indicate  negligible  turning  shear 
for  Case  1  and  moderate  turning  shear  for  Case  2.  Essentially  then,  we  classify  Case  1  as 
a  pure  speed  shear  case,  while  we  classify  Case  2  as  a  speed  and  turning  shear  case.  Case 
1,  whose  cellular  structure  is  aligned  nearly  normal  to  the  mean  wind,  is  fully  consistent 
with  Brown  (1980);  while  Case  2,  whose  rolls  are  aligned  14°  to  the  right  of  the 
background  wind,  is  predicted  to  be  40°  to  the  left  of  the  mean  wind.  However,  Brown 
(1980)  notes  that  in  cases  with  both  speed  and  turning  shear  “rows  may  occur  at  various 
angles,  often  depending  on  the  layer  stratification,  such  that  positive  buoyancy  acts  to 
align  the  rows  with  the  mean  velocity.”  We  are  encouraged  by  these  consistencies. 

4.3.  Case  Study  Conclusions 

The  results  obtained  in  this  study  of  the  Sonde  1  case  fi-om  Hi-Res  2  clearly  show  that 
the  model  is  capable  of  capturing  the  spatial  organization  of  roll  and  cell  circulations  in  the 
MABL.  The  addition  of  an  inversion  layer  appears  to  more  realistically  cap  the  secondary 
circulations  (Laufersweiler  and  Shirer  1995)  than  did  the  rigid  lid  used  in  earlier  runs  and 
in  other  studies  (e.g.,  Lambert  1995).  In  the  presence  of  a  quasi-equilibrated  solution  and 
weak  to  moderate  forcing,  we  can  now  reasonably  expect  the  inversion  to  consistently  cap 
boundary  layer  secondary  circulations.  We  are  also  pleased  that  the  vertical  profiles  of 
heat  and  momentum  fluxes  respond  well  to  the  addition  of  an  inversion  layer,  especially 
matching  the  results  of  Laufersweiler  and  Shirer  (1995)  for  Case  2.  The  appearance  of 
narrow  updrafts  and  broad  downdrafts  in  the  Case  2  vertical  velocity  cross  seaion  is 
another  encouraging  aspect  of  our  findings.  Most  notable,  however,  is  the  ability  of  the 
model  to  reproduce  realistic  stress  patterns  at  the  lower  boundary  and  to  select 
appropriate  orientation  angles  relative  to  the  mean  background  wind.  In  developing  their 


71 


algorithm  to  extract  wind  data  from  SAR  imagery,  Wackerman  et  al.  (1996)  made  the 
assumption  that  stress  patterns  are  inherently  aligned  with  the  direction  of  the  mean 
background  wind  “to  within  a  180®  ambiguity.”  In  most  cases  during  our  study-about  95 
percent  of  the  time—we  found  this  to  be  true  and  therefore  concur  that  they  have  made  a 
valid  assumption.  However,  we  must  always  remember  that  it  is  possible  for  rolls  and 
cells  to  align  normal  to  the  direction  of  the  mean  wind,  especially  in  the  presence  of  an 
inflection  point  in  the  wind  profile;  our  Case  1  results  lucidly  convey  this  result. 

Overall,  the  model  performed  well  in  replicating  the  atmospheric  conditions  of  Hi-Res 
2.  During  the  course  of  our  research,  however,  we  noted  several  areas  that  should  be 
addressed  in  the  near  future: 

1)  After  implementing  the  inversion  subroutine  and  a  numerically  accurate  temporal 
finite  differencing  scheme,  we  did  not  find  attainment  of  a  fully  equilibrated  model  solution 
to  be  possible  within  a  reasonable  amount  of  time.  The  addition  of  a  surface  layer, 
nonlinear  diffusion  and  a  scale  factor  5/ have  helped  tremendously,  but  we  feel  that  further 
numerical  techniques  are  required  to  increase  the  speed  of  integration,  without 
compromising  the  accuracy.  Another  option  would  be  to  find  better  seed  initial 
conditions— using  the  last  vector  _y  from  Case  2  (at  150,000  s)  would  be  a  good  starting 
point  for  all  future  case  studies.  Should  these  proposals  not  resolve  the  problem,  we 
recommend  that  more  vertical  wavenumbers  be  added  to  the  model  to  improve  the 
resolution  given  by  the  spectral  truncation.  The  addition  of  more  horizontal 
wavenumbers,  which  would  allow  the  representation  of  cell  broadening  (Chang  and  Shirer 
1984),  is  warranted  as  well. 

2)  The  subgrid  parameterization  of  fluxes  should  be  addressed  so  that  qualitatively 
correct  heat  and  momentum  flux  calculations  can  be  obtained  down  to  the  lower  boundary 
height  hiB.  This  might  be  accomplished  by  inserting  subgrid  heat  flux  terms  (produced 
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only  by  the  horizontally-independent  modification  coefficients  )  into  the  roll  potential 
temperature  NDS  equations  for  (2.39).  Such  a  parameterization  would  be  an 
application  of  the  eddy  form 


-A. 
a""  ^ 


(w"0") 


(4.6) 


where  the  double  primes  denote  sub-model  values  where  the  number  of  wavenumbers 
exceeds  four  and  is  the  height-only-dependent  portion  of  the  0*  expansion  in  (2.39). 

3)  Study  the  linear  portion  of  the  model  more  fully  to  determine  the  sensitivity  of  the 
preferred  ax  and  Oy  values  to  the  wind  profile,  inversion  strength  and  other  parameters. 

4)  Experiment  further  with  the  scale  factor  Sf,  nonlinear  dissipation  rate  d  and  the 
surface  layer  depth  zsl  to  better  understand  the  effects  they  (individually  and  collectively) 
have  on  the  model  solutions.  An  emphasis  should  be  placed  on  the  physical  significance  of 
these  numerical  adjustments. 

5)  Create  an  algorithm  (to  improve  Woodcock’s  [1975]  diagram)  to  assist  in 
predicting  the  expected  structural  type  of  a  given  regime.  We  need  to  define  how  the 
winds  at  the  top  of  the  dommn  or  how  the  inversion  strength  should  best  be  incorporated 
to  make  improved  planform  predictions,  as  well  as  assess  the  plausibility  of  converting  the 
diagram  to  one  of  probability  type,  on  which  lines  of  equal  probability  of  observing  a 
planform  dimension  are  given.  For  instance,  this  would  help  us  determine  whether  or  not 
our  quasi-two-dimensional  Case  2  is  likely  to  occur  in  the  MABL  with  the  given 
parameter  values,  even  though  it  lies  in  the  cellular  regime. 
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6)  Develop  a  three-dimensional  visualization  algorithm  that  is  incorporated  into  all 
future  model  development. 

7)  Study  further  the  presence  of  vertical  velocity  and  momentum  flux  maxima  in  the 
inversion  layer  for  Case  1. 

8)  Conduct  case  studies  using  Hi-Res  2  Sonde  2  data,  listed  in  Table  4.1,  as  well  as 
other  cases  from  subsequent  field  projects. 

Significant  alterations  have  been  made  to  the  model  during  the  course  of  this  study, 
and  we  are  pleased  that  model  solutions  have  improved  greatly  as  a  result.  With  the 
implementation  of  the  recommendations  noted  above,  the  model  should  become  a  reliable 
tool  to  help  ascertain  the  impact  that  BLCE  have  on  stress  at  the  sea  surface  by  defining 
typical  wind  profiles  and  profile  modifications,  and  by  giving  structural  types  of  stress 
patterns. 
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